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I, Vickie Lien Singleton, declare the following: 

1. I hold a Bachelor of Science degree in Mechanical Engineering from West Virginia 
University, a Master of Science in Civil Engineering, specializing in airlift aerator modelling, 
from Virginia Polytechnic Institute and State University, and am currently a fourth year Ph.D. 
candidate at Virginia Polytechnic Institute and State University studying bubble plume dynamics 
and water quality modelling. A copy of my resume is attached as Exhibit A. 

2. I have read and am familiar with the above-referenced application. I have also read and 
am familiar with the Office Action mailed May 4, 2007 ("Office Action") pertaining to this 
application. 


In re application of: Ronald D. BLUM, ET 
AL. 

Application No.: 09/994,860 


Filed: November 28, 2001 

For: METHOD AND APPARATUS FOR 
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3. I have been retained by the assignee of the above-identified application to assist in 
responding to the Office Action. I have no financial interest in the assignee or the outcome of 
this patent application, including whether it issues as a patent or not. I am being compensated 
for the time spent on this matter at the rate of $200/hr, plus reasonable expenses. 

4. It is my understanding that the claims currently under examination, which relate to, inter 

alia, methodologies for reducing the intensity of a hurricane, were rejected as allegedly being 

wholly inoperative, lacking credible utility, and not enabled. Specifically, the Examiner asserts: 

.. .applicant admits in his arguments, "submersibles of the kind 
required for this application do not presently exist." It seems that 
applicant wishes that someone will come along and develop the 
technology required to make the required submersibles, thereby 
enabling the present invention. Therefore, it is impossible for one 
of ordinary skill in the art at his time to make and or use this 
invention. Office Action at page 6. 

5. This Declaration provides support for an initial order-of-magnitude estimate of the gas 
flow rate and the quantity of gas required to induce an adequate upwelling flow rate to lower the 
temperature of the upper sea surface, using the amount of direction provided for in the 
specification in combination with the knowledge of one skilled in the art at the time of filing the 
application, such as the use of two existing bubble-plume models, without any undue or 
unreasonable experimentation. In particular, the specification is replete with the detailed 
description of various gas sources (see, e.g., para 0029), submersible designs (see, e.g., para 
0034), inception strategies (see, e.g., 0026-27, 49-64), and the exemplary calculation for 
upwelling volume required for successful surface layer water cooling (see, e.g., 0049-64). 

6. The specification is directed to reducing the intensity of hurricanes at sea by upwelling 
deep-water in the ocean. To obtain an initial order-of-magnitude estimate of the gas flow rate 
required to induce an adequate upwelling flow rate, two existing bubble-plume models were 
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applied. Each model is based on the discrete-bubble method developed by Wiiest, A. ET al., 
Water Resour. Res. 1992, 28, 3235-3250, and accounts for ambient density stratification and 
mass transfer between bubbles and water. The circular bubble-plume model was published in 
1992 and is for a circular or round gas diffuser. The linear bubble-plume model was first 
presented in 2001 (McGinnis, D. et al., Asian Waterqual 2001: First IWA Asia-Pacific 
Regional Conference: Fukuoka, Japan, 12-15 Sept. 2001) and was later refined and validated in 
2007 (Singleton, V. et al., Water Resour. Res. 2007, 43, W02405). This model is for a long 
and narrow diffuser. Both diffuser configurations were modeled to compare the effects of 
geometry on upwelled or induced water flow rate. Based on the design example in the patent 
application, it was assumed that the diffusers would be located at about 300 m depth when 
bubbling. The patent application states that carbon dioxide is the preferred gas, so model 
calculations focus on this compound. Because they are used for lake oxygenation studies, the 
solubility (Weiss, R. ET AL., Mar. Chem. 1974, 2, 203-215) and mass transfer coefficients (Clift, 
R., et al., Bubble, Drops, and Particles New York, NY, 1978) of the bubble-plume models 
were modified for carbon dioxide, and the effect of dissolved carbon dioxide on water density 
was included (Weiss ET AL, supra). The Wiiest reference is attached as Exhibit B, the McGinnis 
reference is attached as Exhibit C, the Singleton reference is attached as Exhibit D, the Weiss 
reference is attached as Exhibit E, and the Clift reference is attached as Exhibit F. 

7. To apply the models, boundary condition profiles of temperature and salinity are needed. 
In an attempt to closely represent the design example of the application, profile data collected off 
the eastern coast of Florida (26.17N,-78.80W) on September 19, 2007 were used 
(www.aoml.noaa.gov/phod/triananes/tmp/dxbtl 190840242.dat). The depth of the 26° C 
isotherm in this recent profile is approximately 81 m, which is slightly more conservative than 
the 70 m depth of the patent application design example. The models also require the 
dimensions of the diffusers and the initial bubble size. The initial bubble size was assumed to be 
10 mm. As a starting point, the length of the linear diffuser was assumed to be proportional to 
the length of an SSBN-726 Ohio-Class submarine, which is about 170 m. Using a linear diffuser 
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length and width of 150 m and 13 m, respectively, produces a diffuser area of 1,950 m 2 . An 
equivalent diameter for the circular diffuser for this cross-sectional area is approximately 50 m. 

8. Using these initial and boundary conditions, each model was run over a range of applied 
gas flow rates (2,000-20,000 Nm 3 /s) to a single diffuser system or unit. Because the intent of the 
apparatus is to cool surface waters, the lowest gas flow rate that upwelled water to and detrained 
at the surface was selected. The gas flow rates that induced plumes to the surface were 10,500 
and 12,500 Nm 3 /s for the linear and circular diffusers, respectively. The depth of the 26°C 
isotherm, the temperature below which hurricane development is hampered, was 81 m. 

Therefore, only plume induced flow rates from 300 to 81 m depth are considered as effective 
upwelling flow rates. For the previously stated gas flow rates, the plume water flow rates at the 
depth of 26°C isotherm are 51,900 and 51,300 m 3 /s for a single linear and circular diffuser, 
respectively. These flow rates represent the upwelling flow rate from deeper water into the 
effective epilimnion or surface layer. Even though the plume continues to rise through the 
epilimnion and entrain ambient water, plume induced flow that occurs within the effective 
epilimnion does not contribute to the upwelling into this volume. It would be expected that the 
artificial upwelling of the deep, cold seawater to the sea surface layer by the bubble-driven 
plume would create an upper ocean layer region of sufficiently lower temperature. 

9. Referring to the patent application, the maneuver before upwelling submersible 
mobilization strategy is estimated to require a total upwelled water flow rate of at least 12.1 
million m 3 /s. Using the estimated water flow rate values of 51,900 and 51,300 m 3 /s for a single 
linear and circular diffuser unit, respectively, a total of 233 linear or 236 circular diffusers may 
be needed. This corresponds to total gas flow rates of 2.45 and 2.95 million Nm /s for the linear 
and circular diffusers, respectively. Per the patent application, the maneuver before upwelling 
strategy requires a total gas volume of 5.29 x 10 10 Nm 3 for the linear diffusers and 6.37 x 10 10 
Nm 3 for the circular diffusers. Assuming that 476 Nm 3 of gas would be liberated per cubic 
meter of liquid, the amount of liquid CO 2 required for linear diffusers would be about 1 x 10 8 m 3 
and about 1.3 x 10 8 m 3 for circular diffusers. 
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10. All statements made herein of my own knowledge are true and that all statements made 
on information and belief are believed to be true; and further, that the statements were made with 
the knowledge that willful false statements and the like so made are punishable by fine or 
imprisonment, or both, under section 1001 of Title 18 of the United States Code, and such willful 
false statements may jeopardize the validity of the application or any patents issuing thereon. 

Dat e !I~~ f -T>7 

Lien Singleton 
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10. Education and Experience: 


January 2003 to Present 

■ Enrolled in Ph.D. program at Virginia Polytechnic Institute and State University 

■ Current GPA: 4.0 

June 2000 to January 2003 

■ Employed as Design Engineer at Black & Veatch Corporation, Louisville, KY 

■ Responsible for coordination of air permitting requirements for a $70 million design- 
build project for wastewater solids processing. Coordinated training for plant staff. 
Performed field review for development of as-built drawings. 

■ Detailed design of influent pumping station and associated yard work. Coordinated 
design support groups. Developed hydraulic profile of new and future treatment 
facilities. Responsible for obtaining required permits and approvals. 

■ Assisted with determination of future treatment costs for a 20-year facility plan for a 
130-mgd water utility. Reviewed client’s specifications for distribution system 
projects. Determined storage and feed capacity of chemical equipment. 

■ Drafted preliminary engineering report for $15 million water distribution system 
expansion. Researched permits required for construction of project. 

May 1998 to June 2000 

■ Employed as Staff Engineer at Black & Veatch Corporation, Greenville, SC 

■ Responsible for shop-drawing review and maintenance of construction 
correspondence databases for a $25 million wastewater plant expansion to 7.5 mgd. 

■ Preparation of scope of work, specifications, and cost estimate to provide assistance 
to achieve improved NPDES permit compliance for in-plant sewer overflows. 

■ Served as project engineer for comprehensive facility plan of a 60-mgd plant that 
evaluated unit process operation and expansion to 90 mgd. Responsible for hydraulic 
evaluation of plant, drafting of technical memoranda, and coordination of personnel. 

August 1995 to May 1998 

■ Attended Virginia Polytechnic Institute and State University 

■ Received M.S. in Environmental Engineering 

■ Final GPA: 4.0 (Summa Cum Laude) 

August 1991 to May 1995 

■ Attended West Virginia University 

■ Received B.S. in Mechanical Engineering 

■ Final GPA: 3.6 (Magna Cum Laude) 

11. Publications: [including under maiden name of Vickie Lien Burris] 

Burris, V. L. and Little, J. C. (1998). Bubble dynamics and oxygen transfer in a 

hypolimnetic aerator, Water Science & Technology, 37 (2) 293-300. 

Burris, V. L., McGinnis, D. F. and Little, J. C. (2002). Predicting oxygen transfer and water 

flow rate in airlift aerators. Water Rese arch, 36, 4605-4615. 



Singleton, V.L. and Little, J.C. (2006). Designing hypolimnetic aeration and oxygenation 
systems - A review, Environmental Science and Technology, 40, 7512-7520. 

Singleton, V.L., Gantzer, P. and Little, J.C. (2007). Linear bubble plume model for 

hypolimnetic oxygenation - Full-scale validation and sensitivity analysis, Water 
Resources Research, 43, W02405. 

12. Other projects 

December 2004 to January 2005 

■ Drafted responses to Virginia Department of Environmental Quality questions 
regarding dissolved oxygen criteria for lakes and reservoirs 

March to May 2005 

■ Lake oxygenation modeling expert for oxygenation pre-pilot design investigation for 
Onondaga Lake, NY 


July 2005 

* Co-author author for technical report “Development of Building Blocks to Prescribe 
Ecological Flows for the Rivanna River Watershed” 

December 2006 

■ Second author for technical report “Review of Oxygenation Technologies with 
Special Reference to Application in the Upper Swan Estuary” 

February 2004 to present 

■ Critical peer review of nine manuscripts submitted by other researchers for journal 
publication 
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Bubble Plume Modeling for Lake Restoration 
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Environmental Physics, Institute for Aquatic Sciences and Water Pollution Control 
Eidgenossische Technische Hochschule, Zurich, Switzerland 
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Dieter M. Imboden 
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A steady bubble plume model is developed to describe a weak air (or oxygen) bubble injection 
system used for the restoration of deep stratified lakes. Since the model is designed for two modes of 
operation, i.e., oxygenation and artificial mixing, gas exchange between water and bubbles has to be 
included. The integral model is based on the entrainment hypothesis and a variable buoyancy flux 
determined by the local plume properties and the ambient water column. Fluxes of eight properties are 
described by nonlinear differential equations which can be numerically integrated. In addition, five 
equations of state are used. The model leaves open two initial conditions, plume radius and plume 
velocity. Model calculations with real lake water profiles demonstrate the range of applicability for 
both modes of operation. The model agrees reasonably well with field data and with laboratory 
experiments conducted by various investigators. 


1. Introduction 

In spite of enormous efforts made to fight lake eutrophi¬ 
cation, in many lakes concentrations of phosphorus and 
other planktonic nutrients are still far above critical values 
IRast and Lee, 1983], As a consequence of high primary 
productivity values, oxygen concentrations in the hypolim- 
nion of eutrophic lakes drop to low values or even to zero. 
So-called internal lake restoration measures have been de¬ 
signed to improve the hypolimnic oxygen levels and to limit 
the recycling of phosphorus from the sediments into the lake 
water [ Imboden , 1985]. Two restoration techniques which 
can be used separately or in combination are (!) artificial 
mixing of the water column during the cold season and (2) 
input of oxygen into the hypolimnion during summer in such 
a way as to preserve the stratification. 

Artificial mixing is technically simpler and also cheaper 
and is thus usually the first method to be evaluated. It is 
designed to add extra mixing energy to the lake at the time 
when the density stratification is weakest in order to allow 
water exchange down to the very bottom, or, if mixing 
occurs naturally, to prolong the period of mixing. Complete 
mixing maximizes the uptake of oxygen by natural gas 
exchange at the water surface. In fact, incomplete mixing 
may not only result from low wind exposure but may also be 
a consequence of lake eutrophication: Mineralization of 
biomass in the hypolimnion often leads to the steady accu¬ 
mulation of dissolved substances in the deep waters and thus 
to a chemically induced permanent stratification [Joller, 
1985]. In some lakes, the hypolimnic oxygen reserves may 
not be large enough to prevent the hypolimnion from becom- 
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ing anoxic during the summer, even if complete saturation 
had been achieved during the winter. In such situations, 
oxygen can be artificially introduced into the hypolimnion 
during the summer. This is commonly combined with artifi¬ 
cial mixing during the winter. 

Artificial mixing is achieved, in most cases, by injecting 
compressed air into the deepest part of the lake. However, 
different technical systems have been designed for the injec¬ 
tion of oxygen into a stratified lake. One frequently used 
system is the Atlas-Copco “Limno" device [Bucksteeg and 
HoUfelder, 1978], This system consists of two concentric 
tubes placed vertically in the water column. Air is pumped 
into the lower end of the inner tube, in which the water rises, 
then sinks back through the outer tube and leaves the system 
through outlets placed at the appropriate depth. 

In this article a different approach is described. The 
system “Tanytarsus,” designed by the two Swiss engineers, 
E. Jungo and U. Schaffner, is presently in operation in 
several deep Swiss lakes [Imboden, 1985; Stadelmann, 1988; 
Stockli and Schmid, 1987], The system uses the same lake 
installation for both the injection of compressed air for 
artificial mixing in winter and the injection of pure oxygen 
into the hypolimnion in summer (Figures 1 a and lf>). Pure 
oxygen is used for hypolimnic oxygenation to prevent the 
accumulation of molecular nitrogen which can be toxic to 
fish [Fast and Hulquist, 1982], 

Both modes of operation are physically similar. The main 
difference between them consists in the choice of the initial 
bubble radius and the selection of either natural air or pure 
oxygen. During the winter the air bubbles have to be large 
enough to guarantee that the buoyancy force acts up to the 
lake surface. In contrast, for efficient oxygenation, the 

a vVO** n KuKKIoc Kaiia te\ Ka c moil onoiiflV> t r\ ^iccr>l\ia form. 
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Fig. la. Schematic diagrams of the application of the bubble injection system (below) and of the coupling of the 
plume model with a water quality lake model (above). The two models should run independently of one another, 
alternating over a short time interval (steady input) to supply one another with their respective required input variables. 
Labeled quantities are those appearing in (24) and (25). 


oxygen-rich surface waters. In addition, the bubble plume 
should be weak enough so that the induced water circulation 
is restricted to the hypolimnion. Otherwise, the induced flux 
of nutrients from the hypolimnion into the epilimnion could 
further stimulate primary production in the lake. 

Diffusors for both modes of operation have been success¬ 
fully designed and employed in several medium-sized Swiss 
lakes (Baldeggersee, Sempachersee, Hallwilersee). In order 
to optimize design and operation of the installations with 
respect to oxygen input and energy consumption it is impor¬ 
tant to gain a better understanding of the behavior of the 
bubble plume in a stratified environment. In this article, 
classical plume theory is extended to include gas-water 
exchange of oxygen and nitrogen as well as bubble expan¬ 
sion due to decreasing pressure during bubble rise. The 
mathematical plume model is presented in sections 2.1—2.6; 
then the relationship of the plume model to an overall lake 
model is discussed briefly in section 2.7. In section 3 the 
model is applied to the two standard cases of operation, 
oxygenation and artificial mixing, respectively, and com¬ 
pared with a field experiment. Results and questions remain- 



Mg. 16. lop view ot one amusor unit, consisting ot moes 
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ing open are discussed in section 4, and conclusions are 
given in section 5. 

2. Development of the Bubble Plume Model 
in Stratified Water 

2.1. Previous Work 

Air bubble plumes, induced by the steady release of air 
into water, have been widely studied both experimentally 
and theoretically. Previous work, which has been mainly 
concentrated on plumes in homogeneous water, has revealed 
various characteristics of the flow: the Gaussian radial 
distribution of plume velocity and bubbles, the similarity of 
the radial profiles at different distances from the source, and 
the entrainment of the surrounding water. Measurements 
have been carried out over a wide range of water depths 
from centimeters [Durst et al., 1986] to meters [Kobus, 1973; 
Ditmars and Cederwall, 1974; Goossens, 1979; Fannelop 
and Sjoen, 1980; Milgram and Van Houten, 1982] and up to 
tens of meters [ Topham , 1975; Milgram, 1983] with airflow 
rates from 4 x 10 -7 N m 3 s" 1 (weak plume [Leitch and 
Baines, 1989]) to 0.66 (V m 3 s _1 (strong plume [Topham, 
1975]). 

Theoretical investigations followed the early work of 
Ditmars and Cederwall [1974], who applied the integral 
theory of a single-phase plume [Morton et al., 1956] to the 
two-phase bubble plume. The theory is based on the hori¬ 
zontally integrated equations of conservation of mass, mo¬ 
mentum and buoyancy and on the entrainment hypothesis, 
which states that the volume of entrained water per unit 
height is proportional to both the local plume velocity and 
the plume circumference. Compressibility of the air and 
bubble slipping, i.e., the differential velocity between rising 
bubbles and the plume water, are two features which are 
usually also included in bubble plume models. The effect of 
turbulent transport within the plume was taken into account 
by Milgram [1983] who introduced the momentum amplifi¬ 
cation factor y (ratio of total momentum flux to the momen- 
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WOEST ET AL.: BUBBLE PLUME MODELING FOR LAKE RESTORATION 


3237 


McDougall [ 1978] generalized Ihe integral model for a strat¬ 
ified ambient fluid. Experiments, carried out in strongly 
stratified [ McDougall , 1978] and weakly stratified [Hussian 
and Narang , I983J water revealed a double plume structure 
in which the bubbles remain in the center part, whereas the 
pure liquid outer part of the plume can spread out at certain 
levels. 

If bubble plumes are applied to destratify lakes and 
reservoirs it is important to optimize the airflow rates 
required to overcome a given stratification. A concept of 
maximum efficiency was demonstrated experimentally for 
linear stratification by Asaeda and Imberger [1988], Patter¬ 
son and Imberger [1989] combined McDougall' s [1978] 
bubble plume model with the one-dimensional dynamic lake 
model DYRESM to evaluate destratification in response to 
air bubble plumes and weather. With a similar goal in mind, 
Zic and Stefan [1990] have even developed models to 
simulate the entire flow field induced by a bubble plume in a 
lake. 

In this section, the integral plume theory will be extended 
to include not only bubble expansion [ Ditmars and Ceder- 
wall, 1974], but also gas exchange (dissolution and stripping) 
for the dynamics of a buoyant plume. Previous investigators 
have not included the latter feature, which is important in 
deep lakes and for weak plumes. The bubble model will also 
include the effects of density stratification due to gradients of 
both vertical temperature and dissolved solids. The model 
can cope with arbitrary ambient profiles of dissolved oxygen 
and nitrogen. 

2.2. Model Features and Assumptions 

Previous investigators have presented the basic equations 
for bubble plumes, and discussed the problem of estimating 
suitable parameters such as the entrainment coefficient, a, 
and the spreading ratio of bubbles to fluid flow, A. The new 
features of the proposed model for lake use are (1) finite 
source diameter; (2) induced vertical water velocity at the 
level of the diffusor as an initial condition; (3) exchange of 
gases (oxygen and nitrogen primarily) between bubbles and 
liquid in the plume; (4) conservation equations for specific 
gases in both dissolved and gas phases; (5) variable gas 
composition of bubbles during rise; and (6) arbitrary ambient 
profiles for temperature, salinity (mass of dissolved solids 
per mass of solution), and dissolved oxygen and nitrogen, 
Other features included in the plume model are variable 
bubble radius due to expansion and dissolution; bubble slip 
velocity and gas mass transfer coefficient as a function of 
bubble radius, and solubility constants for oxygen and 
nitrogen as a function of temperature. 

The principal assumptions used in developing the model 
for a steady rising bubble plume in uniform or stratified 
environment are as follows: 

1. The variation of the effective mass density, Ap, is 
important only in the gravity terms and not for mass or 
momentum fluxes (Boussinesq assumption, as in all plume 
theories). This is also equivalent to assuming relative gas 
volume concentration V g « 1. 

2. The mass density of gas is neglected compared to that 
of water in the momentum equation. 

3. Uniform distributions (“top hat”) across the plume 

arf* assumed fnr water velneitv temneratnre dissolved 

solids, dissolved gases and bubble velocity and concentra- 
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Fig. 2. Diagram defining the bubble plume geometry, showing 
circular “top hat” distributions. The top hat radius of the bubble 
core, A b, is smaller than the plume radius b. The latter is assumed 
to apply for all properties of dissolved species, including tempera¬ 
ture (73 and plume velocity (m>). Variables are defined in Tables I, 
3 a and 3 b. 


tion (as well as for individual gases) (Figure 2). This is 
simpler than trying to use Gaussian profiles, and will yield all 
the correct scaling relationships; also for large sources 
whose size is a significant fraction of the depth, we avoid the 
need for special treatment of a zone of flow establishment, 
during which top hat profiles at the source would have to 
change to Gaussian profiles. Morton [1959] and Fannelop 
and Sjoen [1980] have shown that the top hat distribution 
captures the main flow features, and that only the constants 
need to be adjusted. 

4. The plume radius b for dissolved species (dissolved 
solids and dissolved nitrogen and oxygen) and temperature is 
assumed to be equal to the radius of the top hat plume 
velocity, whereas that for the bubbles is only Kb, where A s 
I. (Although it is expected that dissolved substances will 
actually spread slightly faster than momentum, such a dis¬ 
tinction is an unnecessary complication in this case where 
the bubbles are the nuuor source of buoyancy.) 

5. Ambient currents (other than the plume entrainment 
velocity) are assumed to be zero. 

6. The plume is fully turbulent. 

7. The turbulent transport of momentum and scalar 
quantities (heat, salinity, gases) is neglected compared to 
advective transport calculated from the mean flow. This is 
equivalent to taking Milgram' s [1983] momentum amplifica¬ 
tion factor equal to 1.0, which appears reasonable from his 
experiments. 

8. The bubble source is assumed to produce bubbles at a 
constant rate uniformly distributed over the circular source 
area of radius Kb°. 

9. Bubble coalescence is neglected, or else the rate of 
change of the number of bubbles is prespecified. In this 
model, the number flux of bubbles, N, is kept constant with 
height and equal to the number flux at the source. 

10. The bubbles are all of uniform size (or are adequately 

rpnrpspntp.H hv a «m;npn*;inn nf pmial hnhhlpct Tf instpart a 

bubble size spectrum were prespecified, it would be possible 
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TABLE I. List of the Eight Basic Plume Variables 


Variable 

Symbol 

Units 

Vertical velocity of the plume water 

w 

m s _l 

Plume radius for velocity and dissolved species 

b 

m 

Temperature of plume water and bubbles 

T 

°c 

Total dissolved solids in plume water* 

sp* 

kg m " 3 

Concentration of dissolved 0? 

c O 

mol m 3 

Concentration of dissolved N 2 

C N 

mol m ~ 3 

Concentration of gaseous 0 2 t 

m 0 

mol m ~ 3 

Concentration of gaseous N 2 t 

m N 

mol m ~ 3 


Variables with index a , i.e., 7„, S a , c 0a , c Nu , describe “ambient” (i.e., lake water) properties. 
*5 denotes salinity (mass of salt per mass of solution). 

t Amount of gas per unit bulk volume of the bubble-water mixture within the bubble plume of radius 
Xb. 


to replace the number flux of bubbles N by an array 
j*v'(r,, • • ■ , r n ) and to calculate the fate of the individually 
sized bubbles independently. 

11. Water properties of the initial plume are those of lake 
water at that depth (where the diffusor is installed). 

12. Gas exchange between water and bubbles for gases 
other than oxygen and nitrogen (e.g., argon, carbon dioxide, 
methane) is neglected. However, other gases could easily be 
included in the model by adding the relevant conservation 
equations. 

2.3. Conservation Equations 

There are eight basic plume variables (Table 1) to be 
calculated by simultaneous integration of eight differential 
equations describing the plume dynamics based on the 
conservation laws for mass, momentum, and heat (Table 2). 
Several additional parameters are introduced into the model 
which can be calculated either from the corresponding 
equations of state (Table 3 a) or from empirical approxima¬ 
tions taken from the literature (Table 36). In the following 
the reader is referred to Tables 1-3 for the definition of 
variables and parameters. 

Continuity for plume water flow. The rate of change of 
the volume flow is given by the entrainment of ambient water 
into the plume: 

— [-rr 6 J w(l - A 2 VJ] = laixbw (m 2 s“ ! ) (1) 

dz 

where z is the vertical coordinate (positive upward) and a, 
the entrainment coefficient for the top hat plume, is V 2 times 
the corresponding value used for a Gaussian plume model, 
a g , i.e., 

a = y 2 a g (dimensionless) ( 2 ) 


as explained by Morton [1959]. Furthermore, the top hat 
radius is 

b=$b g (m) (3) 

where b g is the plume width in the Gaussian model (w = w m 
exp ( -r 2 /b g )). The term X 2 V g is a correction for the 
volume occupied by the gas bubbles. Since V g < 10 “ 2 for 
bubble plumes of interest for lake restoration, and following 
the Boussinesq assumption, we may rewrite the equation 
without the V g term: 

d , . 

— ivb 1 w) = lernbw (m 2 s ) (4) 

dz 

Momentum flux. We may again neglect the gas volume 
and gas momentum in calculating the momentum flux M = 
7 rb 2 w 2 of the plume. The rate of change of M with height z 
is equal to the buoyant force per unit height of the plume: 

— ( irh 2 w 2 ) = ——— gnb 2 X 2 
dz P p 

+ — —— g-nb 2 (\ - A 2 ) (m 3 s~ 2 ) (5) 

Pp 

( p a and p K are densities of ambient and plume water, 
respectively, and p p is the density of the plume bubble-water 
mixture). On the right-hand side the first term gives the 
buoyant force for the area of the bubble distribution top hat, 
nX 2 b 2 . The second term is the (negative) buoyancy for the 
fluid inside the annulus (Figure 2) between the bubble top hat 
(radius Xb) and the plume velocity top hat (radius b). When 
there is no ambient fluid stratification ( p a = p w ) or when A 
■> 1 , the second term vanishes. 

Heat or temperature flux. Assuming that the coefficient 


TABLE 2. Definition of the Eight Integrated Flux Variables 


Variable 


Formula 


Units 

Volume flux of plume water 

M * 

= Ttb^w 


m 3 s -1 

Momentum flux 

M 

= 776 2 w 2 


m 4 s 2 

Temperature flux (relative to T = 0°C) 

Ft 

= pT 


°C m 3 s 1 

Flux of total dissolved solids 

Fs 

= pSp w 


kg s 1 


D. 

= uc.\i = 0. N 

/ = O, N 

mol s | 

Flux of gaseous 0 2 , N 2 

Ot 

= ttd l A £ {w + H*/,) m it 

mol s 1 
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Parameter 


Partial pressures wilhin 
bubble at depth z 
Total pressure 


Density of bubble-water 
mixture in plume 
Gas volume* 

Ideal gas law 


Van der Waals gas law 
Bubble radius 


TABLE 3 a. Equations of Slate 


Equation 


Pi = + rn N )]p = (G,/(G 0 + G N )]p 

for i = O, N (bar) 

P = Ps + a Si'Pa9 dz = Ps + a p„ gU s ~ z) (bar) 


Water density at atmospheric 
pressure 


Pa = ff 7 U ) + /,(5 U ) for ambient lake water. 
p„. = /r(T) + US) for plume water 


P P = (1 - Vj) p,. 

V g = [(m 0 + m N )/p]/?t 

Implicit equation for V g : [V 9 l[m 0 + m N ) 
• Ip + a A m O + m N YIV g ] » 
r = (3V, A 2 b 2 {w + w h )IA rr/V) 1/3 


Comments 


Total pressure p is given by the hydrostatic 
approximation 

Pi is atmospheric pressure at lake surface 
(bar); z is vertical distance above diffusor 
(m); z, is height of lake surface above 
diffusor (m); g is acceleration of gravity, 
equal to 9.81 m s -2 ; p„ is mean density of 
ambient water; and a is a pressure 
conversion factor (pascals to bars), equal 
to 10 " 5 bar kg -1 m s 2 . 

f r (T) = 999.868 kg m" 3 + 10 " 3 \a,T + 
a : T 2 + a,r’]; o, = 65.185 kg m -J 
"c" 1 ; a-> = -8.4878 kg m" 3 °c~ 2 ; a j = 
0.05607 kg m " 3 °C" 3 ;/,(S) = yS, where 
either y = 0.802 kg m" 3 (%a)~ 1 if 5 is 
salinity expressed as dissolved salt content 
in per mil or y = 0.705 x 10 kg m" 3 
(pS/cm) ~ 1 if 5 is expressed in terms of 
the electrical conductivity measured at 
20°C (modified after Biihrer and Ambuhl 
[1975)). 


R = 8.314 J mol" 1 K 1 (gas constant) and 
t is absolute temperature of plume water 
(K) 

a„ = 1.4 x 10 * bar m fi mol ; b v = 32 
x 10 " 6 m 3 mol" 1 

/V is the number flux of bubbles released at 
diffusor per unit time 


*V, is the gas volume per total volume of the bubble-water mixture in the inner core of the plume (dimensionless). 


of specific heat per unit water volume (p„.k) is constant, and 
neglecting the heat content of the gas volume, the conserva¬ 
tion equation for the heat flux ( nb 2 wp w KT) can be reduced 
to the corresponding expression for the temperature flux F T 
= trb 2 wT : 

— (irb 2 wT) = 2airbwT a (°C m 2 s"') (6) 

dz 

where T a is the ambient water temperature given by mea¬ 
sured temperature profiles or by the lake model. 

Flux of total dissolved solids. For water temperatures 
close to 4°C, the concentration of dissolved solids, ex¬ 


pressed by the salinity S, plays a role in establishing the 
density stratification. Hence it is necessary to include the 
mass balance for dissolved solids: 


— (rrh^wp W S) = 2a itbwp a S„ (kg m 's ') (7) 

dz 

where S a is the ambient salinity. 

Fluxes of dissolved oxygen and nitrogen. There are two 
gases of interest, molecular oxygen (0 2 ) and nitrogen (N 2 ), 
and two phases, dissolved and gaseous, for each. For each 
species the gas transfer term appears as a gain (or loss) in the 
dissolved gas conservation equation and as a loss (or gain) in 


Parameter 


Bubble slip velocity (Figure 5) 


Gas transfer coefficient for molecular 
oxygen and nitrogen (Figure 3) 
Molecular oxygen solubility constant 
(Figure 4) 

Molecular nitrogen solubility 
constant (Figure 4) 

Constants 

Cniiuimiivui vublTIviviH 

Plume radius ratio 


TABLE 36. Parameter Approximations 


Approximation 


w h (r) = vv,(r/r,) 1 357 r/r, < 7.0 x 10"* 

H'i(r) = h -2 7.0 x 10"* < r/r* < 5.1 x I0" 3 

H’j,(r) = w 3 (r/r,) 0 ' 547 r/r, >5.1 x I0 -3 
with W| = 4474 m s _l ; w; = 0.23 ms trj = 4.202 m s 
r, = 1 m 


0,(r) = 

0.6 m s 

1 (r/r.) 

(r/r*) < 6.67 x 10 i = 

O, N 

ms" 1 


ft(r> = 

4 x 10"' 

'ms 1 

(r/r,) > 6.67 x 

!0" 4 . t 

= O, N 

ms" 1 


K 0 iT) 

= K oi + 

k 02 t + 

K m T 2 with K p , = 

2.125 mol m“* 

mol m -3 

bar 

bar 

’> K 02 = 

-0.05021 

mol m 3 bar ' “C 

';k os 

= 5.77 x 



10" 4 

mol m " 3 

bar" 1 °C 

-2 





/Cn(T) 

= ^ni + 

K m? + 

K N3^~ . w > l h ^N1 = 

= 1.042 mol m 

mol m" 3 

bar 

bar 

Km 7 

-0.0245 

mol m bar 1 °C" 

: *n.i 

= 3.171 x 



I0" 4 

mol m 3 

bar" 1 °C 

-2 






u 0.11 

A = 0.8 





3240 


Wuest et al.: Bubble Plume Modeling for Lake Restoration 



Bubbla Rtdlu* [mm] 


Fig. 3. Mass transfer coefficients/3 0 for oxygen measured in tap 
water as a function of bubble radius, based on data (solid circles) 
compiled by Motarjemi and Jameson [1978). In this model the same 
mass transfer coefficients for oxygen and nitrogen were used, since 
the two molecules have nearly identical molecular diffusivities. The 
solid line shows the parametrization /3 0 (r) = 0 N (r) used in this 
paper. 


the gas phase conservation equation. The gas flux through 
bubble surfaces is given by 

Pi(c si - a) (mol m~ 2 s" 1 ) (8) 

where /3, is the mass transfer coefficient (Figure 3) for the gas 
species i (0 2 or N 2 ), c ; is the in situ dissolved gas concen¬ 
tration, and c si is the saturation concentration determined 
by Henry’s law: 

c s , = KiPi (molm~ 3 ) (9) 

Here p t is the in situ partial pressure of the gas phase of 
species i and K t is the solubility constant (Figure 4). 

In order to calculate the total gas transfer per unit height of 
the bubble plume we use the total number flux of bubbles, 
N, a quantity which, according to the list of assumptions 
(see section 2.2, ninth assumption) is determined by the 
initial conditions (initial bubble radius r°, total gas input at 
diffusor and diffusor depth; Tables 3o and 4) and then 
remains constant within the plume. The number of bubbles 
per unit plume height is then simply N/(w + w b ), where (w 
+ w b ) is the total bubble velocity (w b is bubble slip velocity; 
Figure 5). The total surface area of bubbles per unit plume 
height is found by multiplying by the area of each bubble, 
i.e., 4irr 2 NI(w + w b ). The total gas transfer rate per unit 
height is then obtained by multiplying this aggregate surface 
area (per unit height) by the local gas flux; 

4irr 2 N . 

- piiKiPj-Cj) (mol s m" ) (10) 

w + w b 

The conservation equation for each dissolved species t, i.e., 
0 2 and N 2 , includes entrainment of ambient lake water 
(concentrations c ia for i = O and N) and gas exchange: 

d 4irr l N 

— (itb 2 WC 0 ) = 2a TTbwCo a +-;- Poi^oPO ~ c o) 

dz *v + w b 

(mol s’ 1 m~') (11) 

d 4nr 2 N 

— (trb 2 wcti) = lairbwc^ +-/3n(^nPn ~ c 'n) 

dz w r w b 

(mol s _l trT 1 ) (12) 


Fluxes of gaseous oxygen and nitrogen. The rate of 
change of advective flux of gaseous 0 2 and N 2 is equal to the 
rate of dissolution or stripping as determined above: 

d . 4 irr 2 N 

— [irb A (w + w b )m 0 ] = -- Po(KoPo ~ c o ) 

dz w + 

(mol s _l m _l ) (13) 

d 7 . 4nr 2 N 

— [irb A (vr + H’(,)m N ] =--- Ph(KnPn ~ c n) 

dz w + w b 

(mol s -1 m _l ) (14) 

(m 0 and m N are the gaseous concentration of the species per 
unit bulk volume of the bubble-water mixture). Note that the 
relevant volume flux for the gaseous components contains 
the total bubble velocity, w + w b , and the reduced plume 
radius, Kb. Upon solving these equations it is found that the 
proportion of 0 2 and N 2 in the rising gas bubbles is chang¬ 
ing. For instance, pure 0 2 bubbles will strip N 2 from the 
water column at the same time that 0 2 is being dissolved. 

The assumptions made in developing these eight equations 
could be generalized, but at the price of adding parameters. 
For instance, one could distinguish five separate top hat radii 
for velocity, bubbles, dissolved gases, temperature, and 
salinity (giving four different A ratios). The refinement is not 
warranted at this time, since information is lacking on how to 
select different A values. Hence, we use only two top hat 
radii: b for velocity, temperature, salinity, and dissolved 
gases, and Kb for bubbles. It is well known from previous 
investigations that A is significantly less than I, and probably 
about 0.8 [ Milgram , 1983]. 

2.4. Variable Transformation 

To simplify the solution we let the basic dependent vari¬ 
ables be the eight integrated fluxes defined in Table 2, 
namely those for water volume, momentum, temperature, 
dissolved solids and the dissolved and gaseous phases of 
molecular 0 2 and N 2 . With these definitions, the eight 
conservation equations ((4>—(7) and (11M14)) may be rewrit¬ 
ten in terms of the new variables: 



Temperature [°C] 

Fig. 4. Solubility constants K, of molecular nitrogen and oxy¬ 
gen as functions of temperature. Dots represent values taken from 

[ I 07A] QexllA linae -I"- 

in this model. 





Bubble Slip Velocity [cnt/s] 
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TABLE 4. Initial Conditions for Bubble Plume Calculation 


Initial Plume Variable 

Symbol 

Method for Determining Initial Value 

Plume radius 

b° 

effective radius of the entire diffusor unit 
(diffusor geometry) 

Vertical plume velocity 

W° 

h>° depends on source and plume properties: 
see section 2.5 

Volume flux of plume water 


p° = *6°V 

Momentum flux 

M° 

M° = Ttb° 2 w° 2 

Temperature 

r 

ambient temperature at diffusor depth 

Temperature flux 

Dissolved species 

Ft 

F t - n°r 

Concentrations 

S°, c 0 , c N 

ambient concentrations at diffusor depth 

Fluxes 

F’s. Do, o; 

F.s = D 0 = P°Co, Dn = p°c N 

Cas fluxes 

G n , Go 

nitrogen and/or oxygen input rale (mode of 
operation) 

Bubble radius* 

r° 

design (material and porosity) of diffusor 

Pressure 

P° 

depth of diffusor (topography and 
installation) and ambient density profile 

Number flux of bubbles! 

N° 

initial gas fluxes, bubble radius and pressure 
at the diffusor depth (Table 3a) 


’It may be possible to vary bubble radius during operation by switching between different installed 
diffusor tubes. For the system “Tanytarsus" a switch between small (oxygenation mode) and large 
(artificial mixing mode) bubbles is possible. 
tThe number flux of bubbles, N( z), is kept constant in the plume: A/(z) = N°. 



0.1 0.2 0.4 1 2 4 10 20 40 


Bubble Radius [mm] 

Fig. 5. Terminal slip velocity of air bubbles, w b , measured in 
tap water (solid circles) and distilled water (open circles) [Haberman 
and Morion. 19541. The solid line represents the approximation used 
in this model. 



dG Q 4nr 2 N ^ ^ £>o 

~A = ~/lY/ 1 7 Po\ KqPq ' 
dz (Mlp.) + w 'b \ M , 


dG N _ 4irr l N ( £> N 

dz (M/p) + w h p 


( 21 ) 

( 22 ) 


This is a set of eight coupled nonlinear first-order differential 
equations which can be integrated numerically. The system 
contains five additional variables ( p p , p w , r, p 0 , p N ) which 
are linked to basic variables by equations of state (Table 3a) 
and thus can be adjusted after every integration step. Fur¬ 
thermore, there are other physicochemical parameters 
which are either kept constant throughout the model calcu¬ 
lation (e.g., relative bubble plume radius A, entrainment 
coefficient a) or can be calculated from empirical approxi¬ 
mations (Table 3 b). The number flux of bubbles N is 
established by the initial gas volume flux and the assumed 
initial bubble radius, and then remains constant within the 
plume. 


2.5. Initial Conditions 

The model calculations depend crucially on both the initial 
plume conditions, which are summarized in Table 4, and the 
ambient lake water profiles. Whereas the ambient profiles 
are given by the state of the lake (influenced however by the 
plume operation in the long term, see section 2.7), the initial 
conditions can be controlled to a certain extent by the design 
and operation of the diffusor units. 

We have to specify the initial values of the eight flux 
variables defined in Table 2. According to the eleventh 
assumption, the initial values of temperature T°, salinity S°, 
dissolved (J 2 concentration c 0 and dissolved N 2 concentra- 
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tion cn are given by the ambient lake water characteristics at 
the depth of the diflfusor. However, to determine the respec¬ 
tive initial fluxes we need to multiply each by the initial 
volume flux p° which is discussed below. The gaseous 0 2 
and N 2 fluxes, Go and G^, are the key initial variables to be 
controlled by the mode of operation of the diflfusor. The 
initial bubble radius, r°, is determined mainly by the design 
of the diflfusor (pore size) and cannot be controlled by 
varying the gas flow rate (unless it is very large). The number 
flux of bubbles, N = (V°, kept constant in the whole plume, 
is determined by total gas flux, Go + G^, initial bubble 
radius, r°, and pressure, p° (Table 3a), and is not an 
independent initial condition. 

The initial inner plume radius (A b°) is determined by the 
apparent size of the bubble source area. Specifying A and h° 
establishes the width of the velocity plume. In practice, a 
diflfusor unit may consist of several diflfusor tubes distributed 
in a starlike pattern of six arms, as shown in Figure 1 b 
[ Stadelmann , 1988]. In such cases we consider the area of 
the entire diflfusor unit as a uniform source of buoyancy and 
consequently set the initial plume area Tr(Kb°)~ equal to the 
circular area covered by the tubes. 

Mathematically, the vertical velocity w° is an additional 
independent initial variable. In reality, however, the initial 
water velocity should be determined by the properties of the 
source and of the ambient water. It has been found useful to 
define a densimetric Froude number (dimensionless) as 

w 

Fr [2A bg{p a - p p )lp„V n 

For a single-phase plume Fr is simply inversely proportional 
to the local Richardson number Ri [Fischer et a!., 1979] 
which has the unique value of 0.557 along the plume, except 
near the source. The constant of proportionality depends on 
the type of plume profiles being used; here Fr = 
(7 t/4 ) . Thus the Froude number of the plume should 

everywhere have a value of 1.69, except near the source. 
However, for bubble plumes there is an additional charac¬ 
teristic velocity, the bubble slip velocity w h . Thus we do not 
expect a unique starting condition to exist. For very small 
bubbles ( w b —* 0), the initial Froude number should tie the 
same as for a single-phase plume. From the model calcula¬ 
tions that follow we find that as z increases, Fr does indeed 
tend to a nearly constant value of about 1.7. For initial 
values of Fr° < 1.5, the plume water velocity w initially 
(immediately above the diflfusor) increases with increasing z 
(reasonable behavior), while for Fr° > 1.5, w first decreases 
in the model (unlikely), Therefore, for the standard case 
(section 3.1) we use the value in between, namely, Fr° = 

1.6, and for the sensitivity analysis (section 3.2) we analyze 
the effect of varying Fr° values. 

2.6. Solution Procedure 

The stationary solution for the flux variables (Table 2) is 
calculated by a three-step iteration loop. In the first step the 
eight differential flux equations ((I5M22)) are integrated 
using the standard first-order Euler method, starting from 
the initial values listed in Table 4. The vertical coordinate 
increment dz is controlled by keeping the maximum relative 
change of the flux variables below a chosen value (l%o). By 
sensitivitv analvsis it was found that the l%o criterion is 
sufficiently small to guarantee that the final values of the 


plume variables (Table 1) are not increment dependent 
(<l%c relative error). After each iteration the equations of 
state (Table 3a) are solved for the partial pressure of 0 2 and 
N 2 (p o, p N ), the water and plume densities (p H ,, p p ), and 
the bubble radius ( r). In the third and last step, the remain¬ 
ing parameters, bubble slip velocity, w b (Figure 5), mass 
transfer coefficients, (Figure 3), and the solubility con¬ 
stants, K 0 and K s (Figure 4), are calculated from the 
parameterization listed in Table 3b. 

The iteration is interrupted and the program stopped when 
the plume velocity approaches zero or when the plume 
reaches the surface. Comparison of the final plume water 
density with the ambient density profile then allows the 
determination of the equilibrium depth, i.e., the depth to 
which the plume water would sink back if no further mixing 
occurred. However, the fallback plume (which we do not 
attempt to model) entrains some additional water on the way 
down; as a result, in reality it does not fall back completely 
to the calculated equilibrium depth. 

2.7. Plume Model in Relation to Lake Model 

The plume model is only one component of the overall 
lake model. In the simplest case, the lake is considered to 
consist of two parts: the small plume zone (one or more 
vertical columns) and the much larger lake zone outside the 
plumes (Figure la). The lake zone may in many cases be 
considered to be horizontally well mixed or uniform in 
properties (such as temperature T or dissolved 0 2 ); i.e., 
T a (x, z, t ) = T a { z, t). 

The plume is assumed to be in quasi-steady state, operat¬ 
ing for some period of time At (of the order of days) with 
essentially fixed ambient environmental profiles determined 
by the lake zone. Then the lake zone profiles are adjusted to 
reflect the results of the plume action over the given time 
period A/. In other words, the ambient profiles change so 
slowly that they may be considered piecewise constant for 
the purpose of the plume calculations. 

The presentation of the lake zone model is not the purpose 
of this paper. However, for just the physical aspects, the 
basic equations for coupling to the plume model are conser¬ 
vation of water mass and conservation of any scalar species 
or pollutant. For conservation of water mass at level z, the 
vertical rate of change of the vertical flow, wA , through the 
lake cross section A (w is far-field vertical upward velocity, 
Figure la), is equal to the plume outflow per unit height 
q p (z) minus the plume entrainment q e (z): 

d(wA) . 

—— = <7 p~d e (m 2 s ') (24) 

dz 

By definition, either one or the other of q p and q e must be 
zero at a particular level z; i.e., the plume is either entraining 
( q p = 0) or detraining ( q e = 0). 

In the far field of the lake the mass balance of any 
conservative substance with concentration c„(z) is given by 

dc a d(c a wA) , , 

A —+-= < 7 pt' p - q e c a (mol s 'm ') (25) 

dt dz 

where c a and c p describe ambient and mean plume concen¬ 
trations, respectively. 

The plume model uses c a (z) (and T Q (z)) as inputs, and 
nroduces o„. a., and c, as outputs (Ficure la) which are 
used to compute new values of u> and c p at given time 
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TABLE 5. Characteristics of Internal Lake Restoration 
Measures in Baldeggersee (Switzerland) 


Parameter 

Value or Description 


Lake Data 

Location 

47°12'N, 8°16'E, (15 km north of 


Lucerne, Switzerland) 

Surface altitude 

462 m above sea level 

Surface area 

5.2 km 2 

Volume 

176 x I0 6 m 3 

Maximum depth 

66 m 

Oxygenation Mode (May-October) 

Operation 

four to six diffusors at depth 


between 42 and 64 m 

Total Oj input 

3-4.5 x 10 3 kg d _1 

0.025-0.037 N m 3 s~ u 

0 2 input per diffusor unit 

-0.005-0.007 N m 3 s' 1 

Artificial Mixing Mode <November-Aprilj 

Operation 

three to four diffusors or nozzles 

at depth between 58 and 64 m 

Total air input 

0.056 N m 3 s' 1 

Air input per diffusor unit 

0.014-0.017 N m 3 s' 1 


*1 N m ! = I m 3 of gas at I bar and 0°C. 


intervals At by finite difference equations representing the 
lake model (equations (24) and (25)). This procedure was 
demonstrated by Imboden [1985]. 

3. The Case of Baldeggersee: Application and 
Sensitivity Evaluation of the Model 

In this section application of the model to Baldeggersee, a 
highly eutrophic lake (Table 5), is described. Among all 


Swiss lakes, Baldeggersee has the longest record of artificial 
oxygen and air input [Imboden, 1985]. Artificial destratifica¬ 
tion with an air bubble plume began in February 1982 
followed by a first experimental oxygenation of the hypolim- 
nion using pure oxygen during the next summer. 

Vertical profiles of 0 2 , temperature and salinity measured 
in July 1983 (oxygenation) and November 1983 (artificial 
mixing) are chosen as the input data to run the bubble plume 
model (Figure 6). Measured N 2 profiles in Baldeggersee 
were approximately constant and close to the surface satu¬ 
ration value in winter. Hence a constant value of c N = 0.71 
mol m -3 (20 g m -3 ), corresponding to N 2 saturation at T = 
4°C at the altitude of Baldeggersee, was taken as the ambient 
N? concentration. It is assumed that the density of the water 
is defined only by the temperature and the concentration of 
dissolved solids, which consist mainly ofCa 2+ and HCOf. 
The influence of dissolved solids on the vertical density 
gradient is negligible in the thermocline of most lakes, 
although in the prerestoration period chemical density gra¬ 
dients in the deep hypolimnion may be important as well 
[Joller, 1985]. The ambient density profiles for July and 
November 1983 shown in Figure 6 d were calculated from 
water temperature and electrical conductivity using the 
modified equations given by Biihrer and Ambiihl [1975] 
(Table 3a). 

3.1. The Standard Cases 

The model calculations are based on two standard situa¬ 
tions as listed in Table 6. Since the initial values of plume 
size, vertical plume velocity and bubble radius are not 
exactly known, reasonable values for the standard cases 



Density o [kg m" 3 ] Oxygen [mol m 3 ] 


Fig. 6. Temperature, salinity, density (a = p - 1000 kg m 3 ), and oxygen profiles, measured in July (solid squares) 

and November tnne.il r.irelest I OUT in Ftalrteooersee These nmfiles serve as innm Hala for the two standard rase model 
calculations. 
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TABLE 6. Model Input Values for One Diffusor Unit, Applied to Standard Case and Sensitivity 

Analysis 


Parameter 


Standard Case 

Sensitivity 

Analysis 

Dififusor depth, m 

z° 

65 

n.v.* 

Initial plume size, m 2 

■nb°- 

20 

0-50 

Entrainment factor 

a 

Oxygenation 

0.11 

0.05-0.2 

O 2 input,t N m 3 s _l 


0.0062 

0 . 001 - 0.1 

Initial Froude number! 

Fr° 

1.6 <w° - 0.125 m/s) 

1.0-2.0 

Initial bubble radius, m 

r° 

0.001 

I0~ 4 to 

io-' 

Temperature, salinity and oxygen profiles 

Artificial Mixing 

July 1983 (Figure 6) 

n.v. 

Air input,! N m 3 /s 


0.014 

n.v. 

Initial Froude number! 

Fr° 

1.6 (vv° — 0.172 m/s) 

0 . 5 - 2.5 

Initial bubble radius, m 

r° 

0.006 

n.v. 

Temperature, salinity and oxygen profiles 


Nov. 1983 (Figure 6) 

n.v. 


’Here n.v. denotes not varied. 

tl N m 3 = 1 m 3 of gas at 1 bar and 0°C. 

tSee section 2.5 for the choice of the initial Froude number, equation (23). 


were selected followed by a sensitivity analysis of the model 
results with respect to these choices. The calculated vertical 
variations of important model variables along the plume are 
shown in Figure 7 for both modes of operation, i.e., pure 
oxygen input during the summer and air input during the 
winter. 

The most important difference between the two cases 
(Table 6) lies in the initial bubble radius, i.e., 1 mm in July 
versus 6 mm in November. During oxygenation (July) the 


surface-to-volume ratio of the bubbles is large. The bubbles 
dissolve quickly and their volume is only 0.1% of the initial 
value at about 40 m depth (Figure 7c). When the bubbles 
disappear, the colder and heavier water from the lake bottom 
loses its buoyancy (Figure 7c), the plume stops (Figure la), 
disperses, sinks part way, and merges with the surrounding 
lake water at the depth of neutral density. The model 
calculation is restricted to the zone where the plume is still 
rising. At the depth of maximum plume rise (DMPR) some 





Fig. 7. Model calculations for bubble plumes in Baldeggersee using standard parameters from Table 6 and ambient 
nrnfiloe Imm Fionrp f, Profiles of imnortant model variables alone the vertical nlnme axis are shown for the two modes 

of operation: oxygenation during summer (July) and artificial mixing with compressed air during winter (November). 
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Entrainment Factor a [•] 


Fig. 8. Sensilivity of plume model calculated for the case of oxygen input into Baldeggersee using the standard 
situation of July 1983 (Table 6). Only one parameter is varied from its standard value (labeled by “STD") for each 
calculation. See text for definition of equilibrium depth. 


bubbles may still exist if the system operates with large 
enough bubbles. Possible further dissolution or secondary 
plume creation [ McDougall , 1978] by the escaping bubbles 
above DMPR is not considered in this model. If small 
bubbles are used, however, they disappear before the DMPR 
is reached. In the momentum overshoot phase of negative 
buoyancy no secondary plumes have to be considered; the 
plume then behaves as an ordinary one [Morton et al., 1956]. 

In contrast, during artificial mixing with air (November), 
gas exchange is slower since the bubbles have a larger initial 
radius. The rate of increase in bubble size due to decom¬ 
pression is similar to the rate of decrease due to dissolution. 
The net effect, however, is a slight increase in the in situ 
bubble radius (Figure 7c). Therefore, the plume retains its 
buoyancy up to 20 m below the lake surface (Figure If). The 
remaining momentum brings the plume to the surface (Fig¬ 
ures la and lb). 

Note that the differences between the summer and winter 
cases are primarily related to the different initial bubble 
sizes; the chemical composition of the bubbles (oxygen or 
air, respectively) is less important for deep lakes since the 
high hydrostatic in situ pressure makes the gas exchange rate 
nearly independent of the ambient dissolved gas concentra¬ 
tion at this depth. 

The larger size of the bubbles is not the only reason why, 
in the artificial mixing case, the plume reaches the surface. 
Higher gas flux and weaker ambient stratification also sup¬ 
port stronger plumes. Higher gas flux leads to a larger 
density difference (i.e., buoyancy), and a higher initial plume 
velocity (momentum flux). The weak stratification of the 
water column provides little negative buoyancy to stop the 


plume. The effects of initial gas flux and bubble size on the 
plume will be investigated in the following sensitivity anal¬ 
ysis. 

3.2. Sensitivity Analysis 

With respect to the application of the model to internal 
lake restoration the most important model characteristics are 
the vertical extension of the plume, the magnitude of the 
induced vertical volume flux of water, and the fraction of 
oxygen transferred to the dissolved phase. It is therefore 
important to investigate the sensitivity of these quantities 
with respect to those input parameters which either can be 
influenced by the design of the diffusor system or are not 
well known from field observations. It is not easy to compare 
different model results for the winter case in a didactically 
clear way, since, depending on the parameters, the modeled 
plumes may or may not reach the surface. Therefore the 
following sensitivity analysis refers only to the summer case, 
which allows easier comparison (Figure 8). For each calcu¬ 
lation, all parameters except the one which is varied corre¬ 
spond to the standard case (Table 6). 

Among all the parameters, initial bubble size and initial 
gas flow have the largest influence on the behavior of the 
bubble plume, especially on the depth of maximum plume 
rise (DMPR). If the initial bubble radius is less than 1 mm, 
the DMPR is about 25-28 m and is nearly independent of 
bubble size (Figure 8 a). For larger bubbles, the DMPR 
becomes shallower and the plume reaches its highest level 
for a radius of about 1 cm. The curve (Figure 8o) can be 
qualitatively explained in terms of the size dependence of the 
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Fig. 9. Maximum plume rise calculated as a function of initial 
bubble size for the standard situation of July 1983 with and without 
gas exchange (i.e., dissolution of bubbles). Plume dynamics is 
strongly dependent on dissolution for bubbles of radius smaller than 
about 5 mm. 


gas exchange and bubble slip velocity alone. Since the gas 
exchange rate (fractional volume changes per unit time) 
depends on the surface-to-volume ratio of the bubbles, the 
oxygen dissolves faster in the lake if the bubbles are small. 
In addition, smaller bubbles have lower slip velocities and 
thus a longer contact time with the water before reaching the 
DMPR. 

For bubbles with a radius of less than about 0.7 mm, the 
exchange rate (fractional volume changes per unit time) is 
approximately independent of bubble size, since the transfer 
velocity is proportional to bubble radius, whereas the sur- 
face-to-volume ratio is inversely proportional to bubble 
radius (Table 3 b and Figure 3). As a consequence, the 
smaller and slower the bubbles are, the more buoyancy is 
generated and the longer the buoyancy lasts. This causes the 
plume to rise higher. Bubbles with a radius greater than 
about 0.7 mm have decreasing gas exchange rates per unit 
volume (since the transfer rate is constant) and constant slip 
velocities up to a radius of 5 mm. The height of the plume 
will consequently increase with bubble size. For bubbles 
with a radius greater than about 10 mm, the height of the 
plume decreases again since dissolution no longer has any 
influence (Figure 9). The bubbles rise so rapidly that the 
buoyancy (proportional to (w + wj,) -1 ) is significantly 
diminished; the bubbles escape the DMPR and start a new 
plume (not described by the model). 

Because pure oxygen is expensive it is of major impor¬ 
tance to operate a lake oxygenation system so as to minimize 
the amount of oxygen lost by escaping bubbles and to deliver 
the oxygen as exactly as possible to the required depths. As 
shown in Figure 8b, for bubbles with a radius of less than 3 
mm, more than 99% of the O z dissolves. For bubbles with a 
radius of 3 mm, vertical water circulation reaches a maxi¬ 
mum while no 0 2 bubbles escape from the hypolimnion. For 
larger bubbles the total transfer rate decreases rapidly to 
reach a value of less than 10% for a radius greater than about 
3 cm. 

The oxygen input rate directly influences the density and 
thus the buoyancy of the plume. As seen in Figure 8c, an 
increase in the 0 2 input rate from 10 “ 3 to 10 1 N m 3 s 1 
causes the DMPR to rise from 45 m to 10 m depth. For large 
0 2 input rates the DMPR is determined mainly by the 

pncitinn nf the Ihermoo.I in e i.ft the 7nne of maximum 

density gradient. If the total oxygen flow rate and the plume 


intrusion depth range are prescribed, the model allows the 
determination of the minimum number of diffusor units 
needed to reoxygenate at the desired depths. For example, if 
the oxygen input per diffusor is reduced to 50% (e.g., by 
doubling the number of diffusors in the lake), the DMPR 
sinks by 8 m (Figure 8c). Likewise, for artificial mixing with 
air bubbles, the model can be used to find the optimal bubble 
size and the minimum airflow rate as a function of the 
number of diffusors, provided that the ambient stratification 
and the turnover rate are given. 

The initial plume size and initial Froude number (or water 
velocity), two parameters which are difficult to estimate 
without additional and complicated field measurements, do 
not have a significant influence on the DMPR (Figures 8 d 
and 8e). The initial plume size cannot vary more than 
between about half the diffusor area (10 m 2 ) and twice this 
area (40 m 2 ). The initial plume velocity is expected to 
correspond approximately to a densimetric Froude number 
of 1.6, which we choose slightly smaller than the apparent 
asymptotic value in the near field of 1.7. The low sensitivity 
of the model calculation with respect to these two parame¬ 
ters means that during oxygenation the properties of the 
plume are determined mainly by the plume’s dynamics and 
not by the initial conditions. 

As mentioned above, the entrainment factor is an empir¬ 
ical constant which has been determined in several labora¬ 
tory experiments [ Milgram , 1983]. A typical value for Gaus¬ 
sian profiles is 0.08, which is equivalent to 0.11 for top hat 
profiles. Since the entrainment of surrounding water decel¬ 
erates the water plume, the DMPR depends strongly on this 
factor (Figure 8/). 

Bubble dissolution has a drastic effect on plume dynamics. 
In Figure 9 the calculated DMPR is shown as a function of 
initial bubble size with and without gas exchange, respec¬ 
tively. Since small bubbles dissolve faster and have lower 
slip velocities than large ones, the effect of dissolution on the 
DMPR is most prominent for small bubbles. For radii larger 
than about 2 cm, dissolution can be ignored. The result 
dramatically shows the importance of gas exchange for 
bubble plume models, if the initial bubble radius is less than 
about 5 mm, in deep water (65 m depth in this example). 

3.3. Artificial Tracer Experiment 

On August 22/23, 1983, in collaboration with the Institute 
of Geography of the University of Bern, a tracer experiment 
was conducted in Baldeggersee. During a period of 10 min, 
1.5 kg of uranin, a fluorescent dye which had previously 
been dissolved in 200 L of water taken from the lake bottom, 
was introduced into the plume at a depth of 57 m, just above 
one of the diffusors, which at this time was releasing 0.0083 
/Vm 3 s~'of0 2 . The tracer was pumped from a platform into 
a hose fixed to the diffusor unit. The temporal evolution of 
the dye cloud was then measured using an in situ fluorometer 
operated from the platform. 

During the initial phase of the experiment, part of the dye 
rose up to a depth of 5 m. The average rising velocity was 
0.16 m s~' below 40 m depth, and about 0.04 ms -1 above 
this depth (Table 7). After 12 hours, most of the uranin was 
found between 10 and 23 m depth (Figure 10); a small 
amount was found below this layer. 

Since the parameter to which the. model is most sensitive 

namely, the initial bubble radius r°, is not very well known. 
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TABLE 7. Urania Experiment and Model Calculation, August 
1983 



Model* 

Field 

Experiment 

Depth of maximum dye rise 


5 

(above plume), m 

Depth of maximum plume 

13-11 


rise, m 

Depth of tracer insertion, m 


10-23 

Equilibrium depth,) m 

30-25 


Plume velocity, m/s 

57-40 m depth 

0.12-0.13 

0.16 

40-14 rn depth 

0.06-0.08 


40-5 m depth 


-0.04 


‘Parameters: initial bubble radius, 2-4 mm; initial plume velocity, 
0.15 m/s (f>° = 1.6); oxygen flow, 0.0083 N m 3 s -1 ; depth of 
diffusor, 57 m; initial cross-sectional area, 20 m 2 ; initial plume 
radius, 2.52 m; ambient profiles, July 1983 (Figure 6). 

tDepth to which the plume water sinks back once it has lost all 
momentum. The model value was calculated by assuming that no 
entrainment of ambient water into the plume occurs once the plume 
has reached its highest level. Therefore, the value given represents 
a maximum value. 


a rigorous lest of the model is not possible. However, model 
results are found to agree reasonably well with observation if 
r° is taken to lie within the range 2-4 mm (Table 7). Based on 
visual inspection of the bubbles leaving the diffusor tubes 
and on the fact that of the oxygen escaped to the 
surface, the actual size of r° was estimated to lie within this 
range, In the lower part of the plume, the observed upward 
velocity of the dye was close to the modeled plume velocity. 
In the upper part, the upward velocity of the dye was lower 
than the modeled plume velocity, and the dye was found to 
rise to a greater height than that predicted by the model. We 
presume that small secondary plumes generated by undis¬ 
solved bubbles remaining after the main plume has stopped 
rising carried part of the dye up to a depth of 5 m, resulting 
in a relatively low apparent upward velocity (Table 7). The 
fact that the tracer was later detected in the lake mainly 
between the modeled DMPR (11-13 m) and the modeled 
equilibrium depth (30-25 m) supports this interpretation. 


Plume water appears to have mixed with the surrounding 
water at different rates during the sinking phase. Detailed 
tests of the model are planned for newer diffusor systems. 

4. Discussion 

In this section the model is discussed and compared with 
the results of other investigations in order to explore the 
range of its validity and to identify open questions. 

1. The scientific literature contains just a few experimen¬ 
tal data sets which are suitable for comparing model calcu¬ 
lations with measurements. The hydrodynamic part of the 
model was tested by comparing model results with data from 
the Bugg Spring experiment of Milgram [1983]. In this 
experiment air was injected at 50 m depth through a nozzle 
of comparatively small area (about 10~ 2 m 3 ) into a nonstrat- 
ified system. The model was run with entrainment factors 
a, h = 0.116 determined from the observed value a,_and 
adapted to the top hat distribution by putting a lh = v2a g . 
Calculated and observed values for volume flux, plume 
radius and plume velocity agreed within 10% everywhere 
within the vertical profile. The model calculations are also 
close to the proportionality between plume water volume 
flux and the square root of the airflow rate observed by 
Leitch and Baines [1989], Since bubble dissolution was not 
important in their experiment, the bubble solution dynamics 
was tested by comparison with the experiments of Motar- 
jemi and Jameson [1978] in which bubbles of different sizes 
were released at 10 m depth. Calculated and observed 
dissolution rates were equal to within a few percent. 

2. In this model we do not distinguish between the zone 
of flow establishment and the zone of the fully developed 
plume. We also neglect the initial adjustment of the bubble 
slip velocity just above the diffusor. We try to reduce this 
insufficiency of the model by choosing appropriate initial 
plume velocities. There are physical arguments for reducing 
the range of the initial Froude number Fr°. As displayed in 
Figure II for the case of an unstratified water column, the 
Froude number (equation (23)) approaches a maximum value 
of Fr •» 1.7 at about 15-25 m above the diffusor for all initial 
Froude numbers less than 1 .7. In fact, any initial value of Fr° 


Distance from Diffusor [m] 



Depth (m] 


200 


ISO 


100 


50 


Diflusor 


50 



Fig. 10. Distribution of uranin in Baldeggersee 12 hours after injection into the bubble plume just above the diffusor 

at S7 m Henth The two diagrams show trar.er eonrentratinn alono two mutually nemenrlirnlar transect, Moct of the 

uranin is dispersed between 10 m and 23 m depth. (Courtesy of Naturaqua, Bern, Switzerland.) 



3248 


WOEST ET AL.: BUBBLE PLUME MODELING FOR LAKE RESTORATION 




-70 1 —■■■ .. 1 -7<K—— r Y— . -1 

0 10 20 30 40 SO 0.5 1.0 1.5 2.0 2.5 

Water Volume Flux [m 3 /#] Froude Number [-] 


Fig. 11. Plume velocity and radius, volume flux and local Froude number Fr as functions of depth in unstratified 
water calculated with standard parameters for artificial mixing (Table 6) for various initial Froude numbers Fr° ( Fr° 
labeled on curves). Note the contraction and acceleration of the plume for Fr° < 1.0 and the asymptotic behavior of 
Fr for Fr a < 1.7. 


> 1.7 is unrealistic because the plume would have to 
undergo a decrease in velocity immediately above the diffu- 
sor (just where bubbles enter the flow), as can be seen in 
Figure 1 la for the case Fr° = 2.5. For Fr° < 1.5 the plume 
is accelerated immediately above the diffusor, as demon¬ 
strated in Figure 1 la for the case Fr° = 1.0. For Fr° 3 1.0 
the plume initially undergoes contraction as it rises (Figure 
II b). This contraction becomes more pronounced as Fr° 
decreases to 0.5. Because a slight contraction can be ob¬ 
served in tank experiments [Muller et al ,, 1987], we regard 
values of 0.75-1.0 as a reasonable lower bound for Fr°. 

The initial Froude number Fr° depends most likely on the 
geometry of the diffusor source, since for the zone of flow 
establishment it is relevant whether the water can flow 
upward through the diffusor area (open source, diffusor 
above the bed) or whether the water has to flow horizontally 
to the rising plume (closed source, diffusor resting on the 
bed). It seems likely that in the latter case the initial Froude 
number Fr° will be smaller. Because the slight contraction 
which was observed in tank experiments occurred for a 
closed source [Muller et al., 1987] and because the lake 
restoration diffusor unit with widely spaced tubes represents 
an open source, we conclude that the initial Froude number 
for the latter lies between 1.0 and 1.7, 

3. We assumed top hat distributions for all plume prop¬ 
erties. Mathematically, the scaling involved here is nearly 
equivalent to that inherent in the horizontally integrated 
model, presented by Ditmars and Cedenvall [1974]. It sim¬ 
plifies the model, but does not restrict its dynamic general¬ 
ity, as long as the profiles (of velocity and buoyancy) remain 
similar at all heights. We introduced A, the ratio of the inner 
nlnme radius tcontaining the hubbies) to the velocitv radius, 
as the only plume geometry parameter. However, as long as 


the water velocities in the bubble core and the outer annulus 
are equal, the model is independent of A, as is evident from 
(16). We retained A in this model, however, in order to obtain 
the correct equations for 0 2 and N 2 if the double plume 
model [ McDougall , 1978] is used. 

4. In the case of strong stratification, single-phase plume 
models do not adequately describe the plume if the bubbles 
are not completely dissolved at the predicted depth of 
maximum plume rise. Bubble plume experiments carried out 
by McDougall [1978] in a stratified tank showed that fluid 
from the outer annulus can leave the plume at levels of 
neutral buoyancy and spread out horizontally, while the 
center of the plume, where all the bubbles are concentrated, 
continues to rise. The present model disregards the double 
plume structure, the detrainment, and the development of 
secondary plumes initiated by the remaining bubbles. Sec¬ 
ondary plumes are created when oxygenation is carried out 
with relatively large bubbles which are not yet dissolved at 
DMPR, or in the case of artificial mixing, when the airflow 
rate is not large enough to bring the entire plume to the 
surface. For such situations, the model presented here loses 
its validity as soon as the water has lost its vertical momen¬ 
tum, However, oxygen input with adequately small bubbles 
which disappear completely (Figures 7c and le) will not 
generate secondary plumes. The price for using the double 
plume model is the addition of a further two entrainment 
factors. 

5. The entrainment coefficient a is not a constant for all 
plumes. A simple parameterization is not obvious. In our 
model for the standard case a tll = a g = 0.11 was used, 
corresponding to a g = 0.08 as given by Fischer et al. [1979]. 
Mileram f 19831 calculated a for huhhle nlumes of various 
scales based on data from several authors and correlated 
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these values with a “bubble Froude number F B " defined by 
local plume properties. Milgram's [1983] values for a, h 
averaged 0.081 and 0.12 for airflow rates of 0.024 and 0.118 
N m 3 s -1 , respectively, for the 50-m-deep Bugg Spring 
experiment. Fannelop and Sjoen [1980] observed values of 
a,/, from 0.10 to 0.14, increasing with airflow rates in the 
range 0.005-0.022 N m 3 s' 1 . 

6. In the conservation equations used here, vertical 
turbulent flux terms (such as w’c') are neglected compared 
to fluxes based on mean properties, presuming that the 
former are small compared to the latter. This assumption is 
equivalent to a momentum amplification factor of y = 1 
[Milgram , 1983]. The analysis of Milgram [1983, Table 1} 
yielded values of y = 1.0-1.1 for typical oxygenation gas 
rates and y = 1.1-1.26 for gas rates appropriate to artificial 
mixing. The corresponding overestimation of the momentum 
gain would consequently be not more than 21%. However, 
Milgram [1983] has shown that for plumes with low bubble 
density, yean increase dramatically. 

7. There are many questions concerning bubble plumes 
still remaining open which will have to be addressed in tank 
and field experiments. Apart from the problems raised in the 
discussion above, we would like to mention the question of 
how bubble size affects bubble dynamics [Muller et al., 
1987], the problem of continuous detrainment and problems 
related to the interaction of the plume with the ambient flow 
field in the lake (cross flow, plume wandering). Some of 
these questions are being investigated in the new bubble- 
plume tank of the Swiss Federal Institute of Technology 
[Muller et al., 1987]. We think that the mathematical model 
presented here provides a useful framework for further 
investigations and helps to relate laboratory experiments to 
field situations when bubble dissolution and expansion are 
important. 

5. Summary and Conclusions 

1. A steady, horizontally averaged plume model was 
developed based on the entrainment hypothesis, variable 
buoyancy determined by ambient profiles and actual plume 
properties, and processes of bubble gas exchange and de¬ 
compression. The model was applied to a stratified lake for 
a diffusor system with a large open source and low gas flow 
rates (of the order of 5 x 10~ 4 (V m 3 s _l m~ 2 ). Model 
calculations were in encouraging agreement with experimen¬ 
tal data. 

2. Gas exchange, decompression and bubble slip veloc¬ 
ity affect both bubble volume and bubble residence time. 
Consequently, buoyancy flux is a complex function of 
height. For small bubbles, gas exchange determines plume 
dynamics to a large extent. 

3. Apart from the initial bubble size, the gas flow rate 
and entrainment coefficient are the parameters to which the 
model is most sensitive. The sensitivity to diffusor size and 
initial velocity turned out to be small, provided the initial 
velocity was chosen below an upper bound. We argue that a 
reasonable range of the initial velocity is given by the 
condition 1 < Fr° < 1.7 for an open source (diffuser situated 
above the bottom), and Fr° < 1 for a closed source (diffusor 
resting on the bottom). 

4. Model calculations clearly define two modes of diffu¬ 
sor operation: oxygenation (or aeration) and artificial mix¬ 
ing. optimal Duonie size depends crucially on the objectives 


and on the radius-dependent parameters mass transfer and 
slip velocity. In order to prevent oxygen loss by escaping 
bubbles and to achieve insertion into the deepest layers for 
oxygenation, bubbles of radius si mm are the best choice 
(here about 0.8 mm). In contrast, artificial mixing is most 
efficient for an initial bubble radius of about I cm. 

5. Given the operational goals, the model helps to opti¬ 
mize system design in order to conduct the internal restora¬ 
tion of a lake as economically as possible. 

6. The model can be used as a tool to design lake 
aeration systems and evaluate sensitivity to operation pa¬ 
rameters. It can also help in transferring laboratory results to 
the field. Finally, it serves to define areas where additional 
research is needed. 

Acknowledgments. We are grateful to A. Muller, A. Gyr and C. 
Hugi for stimulating discussions. Naturaqua Bern for providing 
Figure 10 and D. Livingstone for improving the English. This work 
was supported by Swiss National Science Foundation grant 20- 
27751.89. 


References 

$ 

Asaeda, T.. and J. Imberger. Structures of bubble plumes in 
stratified environments. Environ. Dyn. Rep. ED-88-250, Univ. of 
West. Aust., Nedlands, 1988. 

Bucksteeg, K., and F. Hollfelder. Versuche zur Bestimmung der 
Leislungsfahigkeit eines Gerates zur Tiefenwasserbeliiftung. 
GWF Gas Wasserfach Wasser Abwasser, 119, 65-72, 1978. 

Biihrer, H., and H. Ambxihl, Die Einleitung von gereinigtem Ab¬ 
wasser in Seen. Schweiz. Z. Hydro!., 37. 347-369, 1975. 

Ditmars .). D.. and K. Cederwall, Analysis of air-bubble plumes, in 
Proceedings, 14th Coastal Engineering Conference, Copen¬ 
hagen, chap. 128. pp. 2209-2226, American Society of Civil 
Engineers. New York, 1974. 

Durst, F., B. Schoenung. K. Selanger, and M. Winter, Bubble- 
driven liquid flows, J. Fluid Mech., 170, S3-82, 1986. 

Fannelop, T. K., and K. Sjoen, Hydrodynamics of underwater 
blowouts, Norw. Marit. Res., no. 4, 17-32. 1980. 

Fast, A. W.. and R. G. Hulquist, Supersaturation of nitrogen gas 
caused by artificial aeration in reservoirs. Tech. Rep. E-82-9, U.S. 
Army Eng. Waterways Exp. Sta., Vicksburg. Miss.. 1982. 

Fischer. H, B., E. J. List, R. C. Y. Koh. J. Imberger. and N. H. 
Brooks, Mixing in Inland and Coastal Waters, Academic Press. 
San Diego, Calif., 1979. 

Goossens, L. H. J., Reservoir Destratification With Bubble Col¬ 
umns, Delft University Press, Delft, Netherlands, 1979. 

Haberman, W. L., and R. K. Morton. An experimental study of 
bubbles moving in liquids, Proc. Am. Soc. Civ. Eng., 80, 379-427, 
1954. 

Hussian, N. A., and B. S. Narang, Simplified analysis of air bubble 
plumes in moderately stratified environments, J. Heat Transfer, 
106, 543-551, 1983. 

Imboden. D. M., Restoration of a Swiss lake by internal measures: 
Can models explain reality? in Lake Pollution and Recovery. 
Proceedings, pp. 91-102, European Water Pollution Control 
Association, Rome. 1985. 

Joller, T., Untersuchung vertikaler Mischungsprozesse mil che- 
misch physikalischen Tracem im Hypolimnion des eutrophen 
Baldeggersees, Dissertation 7830, 94 pp.. Eidgenoss. Tech. 
Hochsch., Zurich, Switzerland, 1985. 

Kobus, H., Bemessungsgrundlagen und Anwendungen fiir 
Luftschleier im Wasserbau, Wasser Abwasser Forsch. Prax., 7, 
168 pp., 1973. 

Leitch, A. M., and W. D. Baines. Liquid volume flux in a weak 
bubble plume, J. Fluid Mech., 205 . 77-98, 1989, 

Marshall, T., Gas Encyclopaedia, p. 1150. Elsevier, New York, 
1976. 

McDougall, T. J., Bubble plumes in stratified environments, J. Fluid 
Mech., 85. 655-672, 1978. 

Milgram, J. H., Mean flow in round bubble plumes. J. Fluid Mech.. 
133, 345-376, 1983. 



3250 


WiiEST et al,: Bubble Plume Modeling for Lake Restoration 


Milgram, J. H., and R. J. Van Houten, Plumes from sub-sea well 
blowouts, in Proceedings of the Third International Conference 
on Behaviour of Off-Shore Structures, edited by C. Chryssosto- 
midis and J. J. Connor, pp. 659-684, Hemisphere, New York, 
1982. 

Morton, B. R,, Forced plumes, J. Fluid Mech., 5, 151-163, 1959. 

Morton, B. R., G. I. Taylor, and J. S. Turner, Turbulent gravita¬ 
tional convection from maintained and instantaneous sources, 
Proc. R. Soc. London, Ser. A, 234, 1-23, 1956. 

Motaijemi, M., and G. J. Jameson, Mass transfer from very small 
bubbles—The optimum bubble size for aeration, Chem. Eng. Sci., 
33, 1415-1423, 1978. 

Muller, A., E. Grass, A. Wiiest, and A. Gyr, Modelling of bubble 
plumes, paper presented at XXII Congress, lnt. Assoc, for 
Hydraul. Res., Lausanne, Switzerland, 1987. 

Patterson, J. C., and J. Imberger, Simulation of bubble plume 
destratification systems in reservoirs, Aquatic Sci., 51, 3-18, 
1989. 

Rast, W., and G. F. Lee, Nutrient loading: Estimates for lakes, J. 
Environ. Eng. N. T,, 109, 502-512, 1983. 


Stadelmann, P., Zustand des Sempachersees, Wasser Energie Luft, 
80, 81-96, 1988. 

Stockli, A., and M. Schmid, Die Sanierung des Hallwilersees, 
Wasser Energie Luft, 79, 143-149, 1987. 

Topham, D. R., Hydrodynamics of an oil well blowout, Beaufort 
Sea Tech. Rep. 33, Inst, of Ocean Sci., Sidney, B. C., 1975. 

Zic, K., and H. G. Stefan, Analysis and simulation of mixing of 
stratified lakes or reservoirs by air bubble plumes, Proj. Rep. 305, 
St. Anthony Falls Hydraul. Lab., Univ. of Minn., Minneapolis, 
1990. 


N. H. Brooks, W. M. Keck Laboratory of Hydraulics and Water 
Resources, California Institute of Technology, Pasadena, CA 91125. 

D. M. Imboden and A. Wiiest, Institute for Aquatic Sciences and 
Water Pollution Control, EAWAG, CH-8600 Dubendorf, Switzer¬ 
land. 


(Received December 4, 1991; 
revised July 6, 1992; 
accepted July 15, 1992.) 



Exhibit C 



Submitted to Asian Waterqual 2001 
First IWA Asia-Pacific Regional Conference 
Fukuoka, Japan, September 2001 


HYPOLIMNETIC OXYGENATION: COUPLING BUBBLE-PLUME 
AND RESERVOIR MODELS 
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ABSTRACT 

Stratification of reservoirs may result in hypolimnetic oxygen depletion with negative consequences for 
water quality in hydropower reservoirs, water supply reservoirs, and cold-water fisheries. Although bubble 
plumes are used to add oxygen to the hypolimnion without significantly disrupting the thermal density 
structure, they nevertheless introduce energy, which causes mixing. The induced mixing changes the 
vertical density gradient, and hence the operation of the plume. To account for this effect, a bubble-plume 
model was coupled in a preliminary fashion with a hydrodynamic reservoir model and then used to predict 
operational data obtained from Spring Hollow Reservoir in Virginia, USA. The coupled model was able to 
accurately predict the extent of hypolimnetic warming. Although oxygen consumption was turned off in 
the reservoir model, the predictions followed the form of the observed oxygen concentration profiles quite 
closely. 
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NOMENCLATURE 

-3 

C dissolved concentration, mol m' . 

2 -1 

E entrainment, m s . 

F d dissolved species flux, mol s' 1 . 

F g gaseous species flux, mol s' 1 . 

F s salinity flux, gs' 1 . 

Ft temperature flux, °C m 3 s' 1 , 

g gravitational acceleration, m s' 2 . 

H Henry’s coefficient, mol m ' 3 bar' 1 . 

Kl mass transfer coefficient, m s' 1 . 

L plume length, m. 

M water momentum, m 4 s' 2 . 

N number flux of bubbles, s" 1 . 

P pressure, bar. 

Q plume flow rate, m 3 s’ 1 , 

r bubble radius, m. 

S salinity, g kg' 1 . 

T temperature, °C. 



v velocity, m s' 1 . 

W plume width, m. 

y gaseous concentration, mol m 3 . 

z depth, m. 

Greek letters 


a entrainment coefficient, 

X spreading coefficient, 

p density, kg m' 3 . 

Subscripts 

a ambient water 

b bubble 

i gas species, oxygen or nitrogen 

p plume water and gas mixture 

w plume water 


INTRODUCTION 

Thermal stratification of lakes and reservoirs can result in substantial hypolimnetic oxygen depletion, 
which may have a negative impact on cold-water fisheries, the drinking-water treatment process, and water 
quality downstream of hydropower reservoirs (Little and McGinnis, 2000). One solution to these problems 
is to install bubble-plume diffusers that replenish hypolimnetic oxygen without destratifying the reservoir 
(Wiiest et al., 1992). While bubble plumes are successful at adding oxygen, the added energy will induce 
some degree of mixing in the hypolimnion, depending on the design and operational characteristics. The 
impact of the induced mixing needs to be considered when these systems are designed. For example, 
partial erosion of the thermocline and warming of the hypolimnion may cause premature destratification of 
the reservoir. Higher temperatures and plume induced mixing may also be responsible for an increase in 
hypolimnetic oxygen demand (Little and McGinnis, 2000). Finally, while plume operation changes the 
thermal structure of the reservoir, the performance of the bubble plume itself depends strongly on the 
vertical density gradient. This complex interaction between plume and reservoir needs to be taken into 
account in the design and operation of bubble plumes. In this paper, a preliminary attempt is made to 
account for this interaction by coupling a model that predicts the performance of a bubble plume with a 
hydrodynamic reservoir model. 

Application of coupled model 

The coupled model was tested using data collected from Spring Hollow Reservoir (SHR) located in 
Roanoke County, Virginia, USA. Spring Hollow is a side-storage reservoir with an existing bubble-plume 
diffuser that uses air as the source of oxygen (Little and McGinnis, 2000). SHR tends to be mesotrophic 
and experiences low dissolved oxygen towards the end of the stratified period. Therefore, it is only 
necessary to operate the diffuser for a period of several weeks in the late summer. Table 1 provides design 
and operational characteristics. Figures 1 and 2 show temperature and oxygen profiles, beginning on Sep 
28 immediately before diffuser start-up. During continuous operation of the diffuser, data were collected 
on Sep 30, Oct 2, Oct 5, and Oct 9. The diffuser was shut down on Oct 13, but unfortunately no data were 
collected on that day. Although the diffuser successfully adds oxygen, it is clear that it also causes gentle 
mixing of the hypolimnion, which in turn causes some warming of the hypolimnetic water. Following shut 
down, additional data were collected on Oct 21 and Oct 28, showing no further warming of the 
hypolimnion, but continued depletion of dissolved oxygen. The coupled model should be able to predict 



the observed rate of oxygen addition as well as the extent of thermocline erosion and hypolimnetic 
warming. 


Table 1. Operating conditions for bubble-plume diffuser in Spring Hollow Reservoir. 


Parameter 

Value 

Maximum depth [m] 

55 

Surface area [10 6 m 2 ] 

0.4 

Total water volume [10 6 m 3 ] 

7.2 

Active diffuser length [m] 

360 

Average diffuser depth [m] 

43 

Air flow rate [Nm 3 Ir 1 ]* 

43 


* 1 Nm 3 denotes 1 m 3 of gas at 1 bar and 0 °C 



Figure 1. Observed temperature profiles 
in Spring Hollow Reservoir. 


Figure 2. Observed oxygen profiles in Spring 
Hollow Reservoir. 


Plume model 

A model for predicting the performance of circular bubble plumes has been developed by Wiiest et al. 
(1992). Based on eight flux equations that are solved simultaneously, the model predicts water flow rate, 
oxygen transfer and plume rise height, given the diffuser depth, applied gas flow rate, initial bubble size, 
and diffuser geometry (Wiiest et al., 1992). A key contribution of this work was the use of a variable 
buoyancy flux to account for gas dissolution as well as decompression. This is especially important when 
pure oxygen is used because the rate of dissolution is very high. The diffuser installed in Spring Hollow 
Reservoir is linear and the eight model equations (Table 2 and 3) were modified to conform to linear, as 
opposed to circular, geometry. The entrainment coefficient, a, was set at 0.08 and the spreading 
coefficient, X, was 0.85, as reported by Fannelop et al. (1991) for linear plumes. The computational 
procedure and other details are the same as for the circular plume model (Wiiest et al., 1992). Although the 
linear plume model has not been fully calibrated, the performance was checked by comparing predictions 
based on a square diffuser with those predicted for a circular diffuser of equal area. Similar results were 
obtained. 



























Table 2. Key model variables. 


Variable 

Formula 

Units 

Entrainment 

E = 2(L + W)av 

m 2 s' 1 

Plume Water Volume Flux 

Q = LWv 

m 3 s' 1 

Momentum Flux 

M = LWv 2 

m 4 s" 2 

Temperature Flux 

Ft = QT p 

°C m 3 s’* 

Dissolved Solids Flux 

F s = QSp w 

gs' 1 

Dissolved O 2 and N 2 Fluxes 

F Di =QCi 

mol s' J 

Gaseous O 2 and N 2 Fluxes 

F Gi =LWX 2 (v + v b ) yi 

mol s' 1 


Table 3. Non-linear differential flux equations. 


Water Volume Flux 

iS = E 

dz 

Momentum Flux 

dM _ p a _Pp gLwx, 2 + Pa ' Pw g LW(l l 2 ) 
dz P p P p 

Temperature Flux 

dF T 

dz= ET - 

Salinity Flux 

7*- = Ep.S. 
dz 

Dissolved Gas Flux 

dFry 47ir 2 N 

—= EC; +—-K L (HjPj -Cj) 

dz v + v b 

Gas Flux 

dFrv 47ir^N 

-r L = —^(HiPi-Ci) 
dz v + v b 


Reservoir model 

CE-QUAL-W2 (W2) is a two-dimensional, laterally averaged, hydrodynamic and water quality model 
(USACE, 1995). W2 includes horizontal, but not vertical, momentum, and accounts for momentum 
transfer from influent streams and the associated shear stress imparted to the surrounding water. W2 was 
applied to SHR for the year 1998 when the aeration system was operated. A temperature calibration was 
performed adjusting only the wind-sheltering coefficient. Based on the model calibration, no significant 
warming is expected to occur in the hypolimnion if the diffuser is not in operation. W2 was then set up to 
run during the aeration period. During diffuser operation there was no inflow into the reservoir and the 
only outflow was seepage loss and water supplied to the water treatment plant. For simplicity, all dissolved 
oxygen (DO) sinks were turned off, with DO treated as a conservative tracer. 

Coupled plume and reservoir model 

The location of the 360-m long diffuser changes in elevation with its deepest point being closest to the dam 
wall. Therefore, as shown in Figure 3, it is represented in W2 as 5 discrete segments in consecutive 
columns. For the purpose of coupling, the plume model is used to compute the flow rate of ambient water 
entrained as the plume rises, and the flow rate, temperature, and oxygen concentration of the water that is 
detrained upon reaching the depth of maximum plume rise. In other words, entrainment removes water of 
known temperature and oxygen concentration from a range of depths, while detrainment returns the sum of 
all the entrained water at the specified discharge location after adding the predicted amount of oxygen. As 
shown in Figure 3, this plume action is simulated within W2 by withdrawing water from the “entrainment” 





cells (marked with X in Figure 3) and discharging it at the “detrainment” cell (marked with O in Figure 3). 
Because of the changing temperature gradient during the course of the simulation, the plume rise height 
also changes. It was therefore necessary in some cases to specify a new detrainment cell during the course 
of the simulation, as indicated in Figure 3. 


Depth (meters) 

Top el Cell 
0.00 
1.90 

3.81 
5,71 

M 7.82 
9.52 
«« 11.49 
iSS 13.93 
15.24 
■ 17.14 
fill 19.05 
20,05 
22.88 
24.78 
28.87 
28.57 
3048 
32 38 
34.29 

38.19 
38.10 
40.00 
41.01 

43.81 
45.72 

47.82 
49.53 
5143 
53.44 

'Cell Width (motors) 
Figure 3. Reservoir grid section for CE-QUAL-W2. 



Although the reservoir is laterally averaged, the entrainment operation is adequately represented using a 
simple withdrawal. However, the entire plume discharge is mixed into the detrainment cell. Because the 
theoretical plume width at the point of detrainment approaches infinity, this instantaneous mixing appears 
reasonable to a first approximation. However, the selection of detrainment cell size might be an important 
consideration when the coupling procedure is refined. 


The “manual” coupling procedure was handled as follows. First, the initial reservoir conditions on Sep. 28 
were used in the plume model to obtain the flow rate of entrained and detrained water. This procedure was 
performed for each of the 5 diffuser segments shown in Figure 3. Then, starting on Sep. 28 at 12h00, and 
using the 5 sets of predicted data, the coupled plume/W2 model was run for 6 hours. The predicted 
ambient temperature and DO profiles obtained at the end of this 6-hour period were then used as input to 
the plume model to generate a new set of predicted plume entrainment and detrainment data for each of the 
5 plume segments. The entire procedure was repeated until Oct. 9 at 12h00. 


RESULTS AND DISCUSSION 

The results of the preliminary coupling procedure are shown in Figures 4 and 5. Excluding the water below 
the diffuser, Figure 4 shows that the coupled model predicts mixing and warming induced by plume 
operation quite accurately. The thermocline is eroded during the entire period of operation, even though 
the plume only partially penetrates the thermocline during the first day of simulation at the 2 highest 
diffuser segments shown in Figure 3. Thereafter, the plume stops between 1 and 3 meters below the 
thermocline. The results suggest that the flow of detrained water causes sufficient mixing near the 
thermocline to result in almost exactly the right amount of warming in the hypolimnion. Factors causing 











mixing include vertical dispersion of horizontal momentum, momentum transfer from the detraining plume 
water, and horizontally generated shear. A more thorough investigation into the dominant mixing 
mechanisms is required. 

Although all oxygen sinks are turned off in the W2 model, Figure 5 shows that the coupled model predicts 
the evolution of hypolimnetic DO fairly well. Once the effects of sediment and other oxygen demands are 
correctly incorporated in W2, the predicted oxygen concentrations will not increase as rapidly. However, 
as currently represented, the plume model under-predicts the rate of oxygen addition because oxygen 
transfer ceases when the plume detrains. Bubbles are still present in the hypolimnion and will continue to 
rise through the water column. In effect, this produces a weaker secondary plume that results in some 
additional mixing and oxygen transfer in the hypolimnion, as well as at higher locations in the reservoir. 
The effects of the secondary plume will be investigated, and if significant relative to the primary plume, 
will be incorporated in the coupled model. 




Temperature (°C) 

Figure 4. Observed versus predicted 
hypolimnetic temperature. 


Dissolved Oxygen (g m" 3 ) 

Figure 5. Observed versus predicted 
hypolimnetic DO. 


FUTURE WORK 

While the preliminary results are very encouraging, the iterative “manual” coupling procedure is time 
consuming. The plume model will therefore be incorporated as a subroutine in W2 so that the computation 
can be done automatically. Further testing of the coupled model will be carried out on new data collected 
at SHR in Virginia as well as at Lake Baldegg in Switzerland. Ultimately, the coupled model could also be 
used to improve the calibration of the plume model. 
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[i] An existing linear bubble plume model was improved, and data collected from a 
full-scale diffuser installed in Spring Hollow Reservoir, Virginia, were used to validate 
the model. The depth of maximum plume rise was simulated well for two of the three 
diffuser tests. Temperature predictions deviated from measured profiles near the maximum 
plume rise height, but predicted dissolved oxygen profiles compared very well with 
observations. A sensitivity analysis was performed. The gas flow rate had the greatest effect 
on predicted plume rise height and induced water flow rate, both of which were directly 
proportional to gas flow rate. Oxygen transfer within the hypolimnion was independent 
of all parameters except initial bubble radius and was inversely proportional for radii 
greater than approximately 1 mm. The results of this work suggest that plume dynamics 
and oxygen transfer can successfully be predicted for linear bubble plumes using the 
discrete-bubble approach. 

Citation: Singleton, V. L„ P. Gantzer, and J. C. Little (2007), Linear bubble plume model for hypolimnetic oxygenation: Full-scale 
validation and sensitivity analysis, Water Resour. Res., 43, W02405, doi: 10.1029/2005WR004836. 


1. Introduction 

[ 2 ] Bubble plumes are used in a variety of industrial and 
environmental applications including mixing in chemical 
reactors, stripping of dissolved gases, containment of spills, 
prevention of ice formation, protection of harbors from 
damaging waves [Fannetep et al., 1991], and destratifica¬ 
tion of lakes and reservoirs [Schladow, 1992]. In addition to 
airlift aerators [Burris et al., 2002] and Speece cones 
[McGinnis and Little, 1998], bubble plumes are commonly 
used for hypolimnetic aeration and oxygenation, which 
preserves stratification of water bodies while adding oxygen 
to the deepest layer. Hypolimnetic anoxia negatively affects 
the drinking water treatment process, cold-water fisheries, 
and water quality downstream of hydropower reservoirs. In 
the United States, releases from hydropower reservoirs 
typically must comply with state water quality criteria for 
minimum dissolved oxygen (DO) concentrations [Peterson 
et al., 2003]. Oxygen depletion may lead to increases in 
hydrogen sulfide, ammonia, and phosphorus and the release 
of reduced iron and manganese from the sediments. Hydro¬ 
gen sulfide, iron, and manganese in drinking water usually 
require additional treatment [Cooke and Carlson, 1989], 
Finally, hypoxia can affect sex differentiation and develop¬ 
ment in fish, resulting in male-dominated populations with 
reduced reproductive success [Shang et al., 2006]. 

[ 3 ] A bubble plume model to predict oxygen transfer 
from linear diffuser systems was presented by McGinnis et 
al. [2001], based on the model for a circular diffuser 
developed earlier by Wiiest et al. [1992], While several 
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models for point-source or circular bubble plumes have 
been proposed [Asaeda and Imberger, 1993; Brevik and 
Kluge, 1999; Ditmars and Cederwall, 1974; Fannelep and 
Sjeen, 1980; Johansen, 2000; Kobus, 1968; McDougall, 
1978; Milgram, 1983; Rayyan and Speece, 1977; Sahoo and 
Luketina, 2003; Schladow, 1992; Wiiest et al., 1992; Zheng 
et al., 2002], less work has been conducted on linear (also 
referred to as line, two-dimensional, or planar) bubble 
plumes. Kobus [1968] developed one of the first analytical 
models for linear bubble plumes, which uses an empirical 
correlation to calculate buoyancy flux. Ditmars and Cederwall 
[1974] presented a model similar to that of Kobus [1968] but 
included bubble slip velocity. Brevik [1977] proposed a 
phenomenological theory for two-dimensional bubble 
plumes comparable with that of Ditmars and Cederwall 
[1974], except that kinetic energy was used to predict 
entrainment. Wilkinson [1979] proposed that full-scale 
linear plumes could be characterized by a Weber number. 
Laureshen and Rowe [1987] presented a model for two- 
dimensional bubble plumes in which plume spreading, 
entrainment, and momentum amplification were assumed 
to be functions of the plume Weber number and empirical 
constants. Fannelep et al. [1991] developed a model for 
linear plumes in shallow water and studied the resulting 
surface currents and recirculation cells. Last, Brevik and 
Kluge [1999] expanded an existing model for linear bubble 
plumes to account for vertical turbulence. Although much 
insight into plume dynamics was gained, none of these 
models for linear or two-dimensional bubble plumes 
accounted for ambient stratification or gas transfer. The 
first linear bubble plume model to include gas transfer was 
presented by McGinnis et al. [2001], who converted the 
circular bubble plume model of Wiiest et al. [1992] to linear 
geometry. The incorporation of gas transfer is critical 
because the rapid dissolution rate of oxygen, and nitrogen 
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when compressed air is used, strongly influences the buoy¬ 
ancy of the plume [Wiiest et al, 1992], Gas transfer is 
especially important in deep water bodies and for weak 
plumes because the increased contact time allows greater 
gas exchange. Last, the prediction of oxygen addition from 
hypolimnetic oxygenation systems is facilitated. Despite 
the usefulness of the linear bubble plume model, it has not 
yet been validated at full scale and over a range of operating 
conditions. 

[ 4 ] Using extensive, high-spatial-resolution conductivity 
and temperature as a function of depth (CTD) transect data 
collected in Spring Hollow Reservoir (SHR), Virginia, 
during diffuser operation in 2003 and 2004, the perfor¬ 
mance of the linear bubble plume model is evaluated. The 
motivation for this work includes verification of model 
performance prior to use for design and investigation of 
critical model parameters through sensitivity analysis. Also, 
the accuracy of model predictions for depth of maximum 
plume rise (DMPR) and induced water flow rate should be 
assessed prior to coupling with lake/reservoir hydrody¬ 
namic and water quality models, such as CE-QUAL-W2 
[McGinnis et al., 2001]. In this paper an improved linear 
bubble plume model is presented, observations and model 
predictions are compared, and results of a sensitivity 
analysis are discussed. 

2. Bubble Plumes in Stratified Water Bodies 

[ 5 ] During hypolimnetic oxygenation with bubble 
plumes, compressed gas is continually supplied to diffusers, 
usually located immediately above the sediments, and is 
allowed to bubble freely. A gas-water plume mixture that is 
less dense than the ambient water is created, which causes 
the mixture to ascend through the hypolimnion. As the 
mixture rises, ambient water is entrained into the plume, and 
the plume width increases. The entrained fluid produces a 
double-plume structure, consisting of an inner core that 
contains the bubble-water mixture surrounded by an outer 
annulus that contains plume water relatively free of bubbles 
[McDougall, 1978]. As the outer annulus entrains stratified 
hypolimnetic water, the plume width increases and the 
density decreases. When the negative buoyancy of the 
entrained fluid exceeds the positive buoyancy imparted by 
the bubbles, the plume detrains water at a rate nearly equal 
to that previously entrained [ Lemckert and Imberger, 1993]. 
At this depth, the velocity of the relatively dense water 
within the plume decreases to zero, and the plume stops 
rising. The detraining plume water then forms an annular 
downward flow immediately outside the outer annulus of 
the upward flowing plume water. The detraining plume 
water entrains ambient water until a depth of neutral buoy¬ 
ancy is reached, where a horizontal intrusion is created into 
the hypolimnion [Asaeda and Imberger , 1993]. The undis¬ 
solved bubbles remaining in the bubble-water mixture 
separate from the inner core flow and continue to rise to 
the surface, repeating the entire process. 

3. Linear Bubble Plume Model 

[6] The linear bubble plume model utilizes the discrete- 
bubble approach, which has also been applied to the airlift 
aerator and Speece cone and was recently reviewed in detail 
by Singleton and Little [2006], The linear bubble plume 


model is composed of horizontally integrated equations 
based on the conservation of mass, momentum, and heat. 
Eight flux equations are solved simultaneously to predict 
water flow rate, plume temperature, oxygen and nitrogen 
transfer and concentration, salinity, and plume rise height, 
given diffuser geometry and depth, applied gas flow rate, 
and initial bubble size (Tables 1 and 2). The model accounts 
for density stratification due to vertical temperature and 
salinity gradients. Entrainment is assumed to be propor¬ 
tional to the local (with respect to depth) plume water 
velocity and perimeter. Bubble size varies as the bubbles 
rise due to expansion and dissolution, and bubble slip 
velocity and gas transfer coefficients are functions of 
bubble radius [Wuest et al., 1992]. Also, Henry’s constants 
for oxygen and nitrogen are functions of temperature 
[Wiiest et al., 1992]. The bubble plume mode! equations 
were originally developed by Wiiest et al. [1992] for 
circular geometry but were modified by McGinnis et al. 
[2001] for the linear geometry of the system installed in 
SHR. The equations that include the spreading coefficient 
(A) were recently refined [ Singleton and Little, 2005] to 
more accurately reflect the geometry of the plume at the 
ends of the linear diffuser (Tables 1 and 2), which is 
approximated in plan view as a long, thin rectangle. 
Additional refinements to the model of McGinnis et al. 
[2001], which are detailed in the following paragraphs, 
include use of a correlation to calculate initial bubble size 
[McGinnis and Little, 2002], correction of the entrainment 
coefficient (a) and A for top-hat profiles, and use of a 
Froude number (Fr) to calculate initial water velocity 
[Fischer et al., 1979; Wiiest et al., 1992], Also, water 
quality profiles from the plume near-field, as opposed to 
the reservoir far-field, were used as boundary conditions 
[McGinnis et al, 2004], 

[ 7 ] The entrainment coefficient and A were set at 0.11 and 
0.93, respectively. These values were derived by Fannelfip 
et al. [1991] by fitting Gaussian profiles to laboratory data 
and were modified for the top-hat profile assumption of the 
model using [Fannelep and Sjoen, 1980] 

oct — s/lac ( 1 ) 



where the subscripts T and G refer to top-hat and Gaussian 
profiles, respectively. For simplicity, top-hat or uniform 
profiles are assumed for water velocity, temperature, 
salinity, dissolved and gaseous constituents, and bubble 
velocity [Wiiest et al, 1992]. Other model assumptions are 
as follows: (1) The linear plume width W for temperature 
and dissolved constituents is equal to the width of the plume 
velocity profile, whereas the bubbles are confined to an 
inner core of width \W (A < 1); (2) ambient currents are 
negligible; (3) the diffuser produces bubbles at a constant 
rate and uniform size that are evenly distributed over the 
cross section of initial width A W„; (4) bubble coalescence is 
neglected; (5) initial water properties of the plume are those 
of ambient water at the diffuser depth; and (6) exchange of 
gases other than oxygen and nitrogen is not considered. 

[8] The model predictions are strongly dependent on the 
initial plume conditions and the plume boundary conditions 
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Table 1. Key Variables of the Linear Bubble Plume Model 3 


Variable 

Formula 

Units 

Entrainment factor 

E = 2(L+ W)a v 

m 2 /s 

Plume water volume flux 

Q = LWv 

m 3 /s 

Momentum flux 

M = LWv 2 

m 4 /s 2 

Temperature flux 

Ft = QT P 

°C m 3 /s 

Dissolved solids flux 

F s = QSp„ 

kg/s 

Dissolved 0 2 and N 2 fluxes Fd = QQ 

mot/s 

Gaseous 0 2 and N 2 fluxes 

F Gi = XW{L - W(\ - 

A)] (v + mol/s 


“Revised from McGinnis el al. [2001], 


of temperature, dissolved oxygen, and salinity. Initial con¬ 
ditions for the bubble plume model were determined as 
detailed by Wiiest et al. [1992], except for the following 
deviations. Initial bubble size is calculated using the corre¬ 
lation developed by McGinnis and Little [2002] for the type 
of linear diffliser installed in SHR: 

rf 3 ,2 = 1-12 + 0.938?. (3) 

The correlation was determined using measured Sauter- 
mean bubble diameter (dip) values of 1.1-2.2 mm 
collected over actual unit gas flow rates q of 0.08- 
0.88 m 2 /h at the diffuser. 

[ 9 ] For the circular bubble plume model, Wiiest et al. 
[1992] proposed that the induced vertical water velocity v 
at the diffuser depth is equivalent to the initial plume 
water velocity. To estimate the initial velocity, Wiiest et al. 
[1992] defined a densimetric Fr and utilized a relationship 
between the local Richardson number ( Ri ) and Fr derived 
by Fischer et al. [1979] for single-phase, round buoyant 
jets discharging vertically. The circular bubble plume 
model was recently validated by McGinnis et al. [2004], 
so a similar procedure was used to determine the initial plume 
velocity for the linear bubble plume. Corresponding equa¬ 
tions for Fr were derived for planar or linear plumes using 
relationships presented by Fischer et al. [1979] to obtain 

Fr = Rn (4) 


rr - --T, w; 

[>~Wg[p a - P p )/Pp]' 

where g is gravitational acceleration, p a is the ambient water 
density, and p p is the bubble plume density. The local Ri for 
planar jets has a constant value of 0.735 at distances from 


the source where the flow is more like a plume [Fischer et 
al, 1979]. Consequently, Fr for a planar or linear plume is 
equal to 1.26, except close to the source. Wiiest et al. [1992] 
assumed that the bubble slip velocity near the source was 
relatively low, so the initial Froude number ( Fr „) should be 
equal to the value for a single-phase plume. Unlike the 
circular plume model, the Fr profiles predicted by the linear 
bubble plume model continually increase with depth for the 
diffuser installed in SHR (not shown). For Fr a < 2.0, the 
plume velocity initially increases with decreasing depth 
immediately above the linear diffuser (Figure la). This 
effect was also predicted by Fannelep and Webber [2003] 
for buoyant plumes rising from areal sources, where a point 
of maximum velocity occurred above the source. Addition¬ 
ally, the plume neck (point of minimum radius or width) is 
always below the point of maximum velocity in the plume 
[Fannelep and Webber, 2003], A neck is not predicted for 
the linear bubble plume when Fr 0 = 2.0 (Figure lb). 
However, necking or contraction will only occur when the 
momentum immediately above the source is relatively low 
and/or the entrainment coefficient is relatively low 
[Fannelep and Webber , 2003], Also, Wiiest et al. [1992] 
reasoned that a plume from an open source (diffuser above 
sediments) may not contract because initially entrained 
water is not obstructed. Therefore Fr„ for an open source 
will likely be higher than that for a closed source (diffuser 
resting on sediments). Because the linear diffuser in SHR 
is an open source (Figure 2), Fr„ for the linear bubble 
plume was assumed to be 2.0, and a sensitivity analysis 
was conducted to determine the effect of varying Fr„. 

[ 10 ] The differential flux equations of the linear bubble 
plume model (Table 2) were solved numerically using the 
fourth-order Runge-Kutta method. Further information on 
the general solution procedure, equations of state, and 
model assumptions is provided by Wiiest et al. [1992] and 
McGinnis et al. [2004], The model calculations are only 
valid over the plume rise height, up to the DMPR. When the 
plume stops rising, a secondary plume may form above as 
bubbles that are not completely dissolved continue to rise 
[Asaeda and Imberger, 1993; McDougall, 1978; Schladow, 
1992]. This phenomenon can occur when a bubble plume is 
released into strong density stratification. 

4, Application to Spring Hollow Reservoir, 
Virginia 

4.1. Field Data Collection 

[ 11 ] To fully evaluate the linear bubble plume model, 
experimental data for boundary conditions, rise height, and 
in-plume constituent profiles are required. Testing was 


Table 2. Nonlinear Differential Flux Equations of the Linear Bubble Plume Model 3 


Flux 

Equation 

Water volume 

dQ/dz = E 

Momentum 

dMtdz = [(p„ - p, v )/ppl?LIF+ [(p,„ - pf,)/p p ]g\W[L - W( 1 - A)] 

Temperature 

dFjhlz = ET„ 

Salinity 

dFJdz = EftA 

Dissolved gas 

dF n !dz = EC„ + [4 -kFNKv + v h )]K L (H,P, - C.) 

Gas 

dF a jdz - —[4nr 2 M(v + v*)] K^Hf, - C,) 


“Revised from McGinnis el al. [2001]. 
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Figure 1. Effect of initial Froude number (Fr a ) on 
predictions of plume water velocity and plume width using 
linear bubble plume model. 

conducted using a full-scale linear diffuser (Figure 2) 
installed in Spring Hollow Reservoir, Virginia (Figure 3). 
Constructed in 1995, SHR is a small monomictic, mesotro- 
phic side-stream reservoir that is generally stratified from 
May to December. The reservoir is managed by the Western 
Virginia Water Authority and serves as one of the principle 
drinking water sources for Roanoke County. The water 
body has a maximum depth of 65 m and a maximum 
surface elevation of 431 m. The approximate surface area 
and volume are 0.54 km 2 and 12.4 x 10 6 m 3 , respectively. 
To prevent anoxia in the hypolimnion and the associated 
deterioration of raw water quality, a linear diffuser equipped 
with fine-bubble porous hoses was installed in 1997 
(Figure 3). The 305-m-long diffuser can be supplied with 
compressed air or pure oxygen at various gas flow rates and 
is located in the deepest portion of the reservoir (368— 
372 m elevation). On the basis of an average surface 
elevation of 430 m, the depth of the diffuser during testing 
ranged from 58 to 62 m along its length. Table 3 provides 
the diffuser operating parameters in 2003 and 2004. 

[ 12 ] Diffuser tests were performed in 2003 using com¬ 
pressed air (21% 0 2 ) supplied at a high gas flow rate 
(45 Nm 3 /h average) during 29 June to 14 July and pure oxygen 
(97% 0 2 ) supplied at a low flow rate (11 Nnr/h average) 
during 14-26 August. A third test was conducted in 
2004 using pure oxygen, but at a higher gas flow rate 
(40 Nm 3 /h maximum) during 22 October to 5 November. In 
2003, the diffuser quickly mixed the rather small hypolim- 
netic volume during both tests, and the water quality con¬ 
ditions on the dates of data collection (2 July and 17 August 
for compressed air and pure oxygen, respectively) were by 
that time relatively homogeneous as a result of plume- 
induced mixing. One of the primary objectives of the 2004 


experiments was to maximize the plume signature in the 
hypolimnion and to increase confidence in the linear 
plume model validation under a different set of boundary 
conditions. In 2004, the data were therefore collected on 
24 October, soon after start of diffuser operation. Addition¬ 
ally, the 2004 test was performed later in the stratified season 
to maximize the ambient DO and temperature gradients in 
the hypolimnion. 

[ 13 ] To establish appropriate boundary conditions for the 
plume model, characterization of the plume near-field 
environment is necessary [McGinnis et al., 2004], Therefore 
the data collected included numerous high-spatial-resolution 
CTD (Sea-Bird model SBE 19plus; 4-Hz sampling rate) 
transects measured almost daily before, during, and after 
diffuser operation. The CTD profiler was also equipped 
with a DO probe (1.4-s response time measured at 20°C). 
Profiles were obtained laterally across the diffuser at 0.5-m 
increments for 0-10 m, 2-m increments for 10-20 m, and 
5-m increments for 20-40 m from the centerline of the 
diffuser in both directions (Figures 3 and 4). (The diffuser 
centerline location is shifted to the left for 2004 (Figures 4e 
and 4f) because the diffuser was repositioned earlier in the 
year. Also, the operational length of the diffuser was 
decreased for that year (Table 3).) 

[ 14 ] The geometry of SHR affects the extent and rate of 
circulation within the hypolimnion induced by the bubble 
plume. Because of the relatively small size of SHR, oper¬ 
ation of the linear diffuser created uniform conditions below 
the thermocline within days after startup. While SHR 
bathymetry influences plume-induced mixing in the hypo¬ 
limnion, the effect on short-term plume operation is negli¬ 
gible because the time that individual bubbles spend in the 
hypolimnion is of the order of minutes. The bubbles and 
resulting plume experience a pseudo steady state with 
respect to ambient conditions. Effects on data collection 
due to the ends of the linear diffuser were assumed to be 
negligible because the lateral profile location was over 
150 m from a diffuser end (Figure 3), and the diffuser is 
designed to release a uniform gas flow along its length. 

4.2. Observations and Model Validation 

[15] Plume rise height, spreading, and constituent profiles 
predicted by the linear bubble plume model were compared 
with experimental observations. Critical model input param¬ 
eters for the three test conditions (2 July 2003, 17 August 
2003, and 23 October 2004) are shown in Table 3. The 
boundary conditions were obtained from averaged near¬ 
field lateral profiles [McGinnis et al., 2004] (±2 m and ±1 m 
from plume centerline for 2003 and 2004, respectively) and 
differed significantly between 2003 and 2004 (Figure 5). 
(The presence of two thermoclines in 2003 is due to the 
pumped reservoir inflow that discharges at 396-m elevation, 
which corresponded to approximately 34-m depth during 
diffuser testing. The lower thermocline delineates the effec¬ 
tive hypolimnion for the oxygenation system.) 

[ 16 ] Measured contours of temperature and DO are 
shown in Figure 4, along with corresponding model pre¬ 
dictions for plume width and the DMPR. The actual plume 
boundaries are not well defined in the contour plots, so 
comparison with predicted plume widths is difficult. The 
lack of distinct plume boundaries was due to the almost well 
mixed conditions in the hypolimnion as a result of diffuser 
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Figure 2. Photograph and schematic of linear bubble 
plume diffuser in Spring Hollow Reservoir, Virginia. 
Courtesy of Mark Mobley, Mobley Engineering, Inc. 


operation, particularly for 2003 (Figures 4a-4d). Also, use 
of compressed air did not produce a strong DO plume 
signature for the July test compared with the August and 
October tests with pure oxygen. The actual plume rise 
height is easier to distinguish, especially in the DO contours 
for August and October. The predicted depths of maximum 
plume rise are 38.2, 46.1, and 48.8 m for July, August, and 
October, respectively. For July and October, the DMPR is 
simulated well by the model (Figures 4a, 4b, 4e, and 4f). 
However, the model appears to underestimate the plume rise 
height for August, when the gas flow rate was compara¬ 
tively low (Table 3). The under-predicted DMPR may have 
been due to an overestimated value for a. In a detailed study 
of round plumes, Milgram [1983] found that a is directly 
proportional to the plume gas holdup or fraction. However, 
the linear bubble-plume model assumes that a is constant 
(Table 3). 

[n] The structure of the plumes is similar to those 
observed by Asaeda and Imberger [1993] for round bubble 
plumes in weak stratification (Figure 4). Depending on the 
gas flow rate and stratification strength, three types of 
horizontal intrusions from the plume were reported. The 
pattern of the DO contours for October (Figure 4f) closely 
resembles type 1, which corresponds to a high gas flow rate 
or weak stratification [Asaeda and Imberger , 1993], This is 
similar to a single plume impinging on a free surface, which 


is analogous to the thermocline in SHR. Although the 
plume for July is not easily discerned from the temperature 
or DO contours, the structure is most likely similar to a type 
1 plume. The July plume appears to detrain primarily at the 
lower thermocline with one strong intrusion, as evidenced 
by accumulation of higher oxygenated water near the top of 
the plume (Figure 4b). The plume for August is best 
classified as type 3, which is for low gas flow rates or 
strong density stratification [Asaeda and Imberger , 1993]. 
Type 3 plumes do not have steady intrusions, but instead are 
characterized by alternating, collapsing eddies that cause the 
plume to meander. 

[is] Referring to the October test (Figure 4f), the higher 
DO concentrations at lower depths adjacent to both sides of 
the plume were likely the result of detrained water that sinks 
past the equilibrium depth due to momentum. The temper¬ 
ature isotherms were also depressed immediately beside the 
plume (Figure 4e). This phenomenon was also observed by 
McGinnis et al. [2004], The DO concentration immediately 
above the predicted DMPR and near the vertical plume 
centerline for October 2004 is higher than the ambient 
concentration (Figure 4f). This could have been caused by 
the formation of a secondary plume above the DMPR 
resulting from incompletely dissolved bubbles. The model 
estimates that the bubble size at the top of the first plume 
was about 5 x 10 -5 m for October 2004, which is relatively 
large compared with July 2003 and August 2003 (Figure 6f). 
These undissolved bubbles could have created a secondary 
plume, which entrained oxygenated water from the detrain- 
ment of the first plume and carried it higher into the water 
column. 

[ 19 ] Vertical profiles of constituents and properties within 
the plumes were also predicted for the three diffuser tests 
(Figure 6). For the July test with compressed air, the higher 
gas flux creates a greater buoyancy flux and a higher initial 
plume water velocity (Figure 6d). Also, the concentration 
driving force for oxygen transfer is lower compared with 
pure oxygen, which decreases the rate of bubble dissolution 
with depth (Figure 6f). These effects result in a higher 
DMPR for July compared with August (Figures 4 and 6). 
The model predictions for 2004 differ from those of 2003 
because of the differing boundary conditions (Figure 5). The 
plume rise height for the October test with pure oxygen is 
less than that for August, even though the gas flux was more 
than tripled (Table 3). The ambient temperature, and hence 
density, stratification was stronger in 2004, which provided 
greater negative buoyancy to decrease plume momentum. 
The stronger ambient density stratification also caused the 
plume velocity to decrease more rapidly with depth in 
2004 despite a higher initial velocity from the diffuser 
(Figure 6d). The lower plume rise height in October resulted 
in lower plume water flow rates than July (Figure 6e), even 
with comparable gas fluxes applied (Table 3). 

[ 20 ] Average temperature, DO, and density profiles 
within the plumes were also measured (Figures 6a-6c). 
The temperature, and consequently density, predictions for 
July and October deviate from the measured profiles where 
the plumes reach the top of the hypolimnion (Figures 6a 
and 6c), or where the rate of plume spreading is greatest 
(Figure 4). The model underpredicts the final plume 
temperature by approximately 0.3° and 0.2°C for July 
and October, respectively. One reason for the discrepancy 
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Figure 3. Bathymetric map of Spring Hollow Reservoir, 
Virginia, showing locations of linear bubble plume diffuser 
and conductivity-temperature-depth (CTD) lateral transects. 


may be that the model assumes a is constant. In addition 
to the dependence on plume gas holdup, Milgram [1983] 
also found that a for a circular plume is directly propor¬ 
tional to the local plume radius. As the plume approaches 
its maximum rise height, the width increases rapidly 


(Figure 4). If linear plume dynamics are similar to circular 
plumes, then a for linear plumes may also increase as the 
plume width increases. Additionally, the boundary profiles 
selected may not accurately reflect the actual ambient 
conditions immediately adjacent to the plume along its 
entire rise height. Measured vertical CTD profiles at the 
estimated plume width were averaged and used for boundary 
conditions. However, the plume width varies greatly with 
depth (Figure 4), so use of vertical profiles at a single lateral 
distance from the diffuser is not appropriate. Also, the 
average plume width was visually estimated from the 
temperature and DO contour plots (Figure 4), but the actual 
plume boundaries are not well defined. 

[ 21 ] The model predicts the plume DO profiles well for 
all three diffuser tests (Figure 6b). For July and August the 
model characterizes the initial increase in DO immediately 
above the diffuser quite accurately. The initial rapid increase 
in DO for August and October is due to the higher oxygen 
saturation concentration at depth with the use of pure 
oxygen. The high hydrostatic pressure causes the oxygen 
transfer rate to be almost independent of the ambient DO 
concentration at the depth of diffuser. By contrast, the initial 
DO increase for the July test with compressed air is more 
modest (Figure 6b). The shape of the predicted DO profile 
for October differs somewhat from the experimental data, 
even though the final values at the DMPR differ by only 
0.3 g/m 3 . The hypolimnion and hence the plume near-field 
were more heterogeneous in 2004 than 2003, which may 
have contributed to the overprediction at lower depths if 
the selected boundary profiles did not accurately represent 
water entrained into the plume. 

[ 22 ] Another source of inaccuracy could be the correla¬ 
tion equation used to calculate initial bubble size (equation 
(3)). This relationship was developed using data collected 
over actual air flow rates per unit length of diffuser of 0.08- 
0.88 m 2 /h [McGinnis and Little, 2002], The actual gas flow 
rates per unit length of diffuser for the July, August, and 
October tests were 0.019, 0.0063, and 0.024 m 2 /h, respec¬ 
tively. Therefore the initial bubble diameters used in the 
model were extrapolated beyond the valid correlation range. 


Table 3. Conditions for Linear Bubble Plume Model Validation and Sensitivity Analysis for Spring Hollow Reservoir, Virginia 2 


Parameter 

2 July 2003 

17 Aug 2003 

23 Oct 2004 

Sensitivity Range 

Oxygen in gas supply, % 

21 

97 

97 

n/a 

Entrainment coefficient f ] 

0.11 

0.11 

0.11 

0.05-0.2 

Spreading coefficient [ ] 

0.93 

0.93 

0.93 

0.5-1.0 

Initial Froude number [ ] 

2.0 

2.0 

2.0 

1.0-3.0 

Operational diffuser length, m 

300 

300 

250 

60-625 

Initial plume width, m 

0.16 

0.16 

0.16 

n/a 

Initial plume area, m 2 

50 

50 

42 

10-100 

Gas flow rate, 1 * Nm 3 /h 

38 

13 

40 

1-400 

Initial gas fluxj m/h 

0.76 

0.26 

0.96 

0.02-9.5 

Initial bubble radius, m 

5.7 x 10“ 4 

5.6 x 10“ 4 

5.7 x 10 4 

10 4 -l0 1 

Diffuser depth, d m 

60 

60 

59 

n/a 

Reservoir maximum depth, m 

64 

65 

64 

n/a 

Reservoir surface area, 10 6 m 2 

0.53 

0.53 

0.53 

n/a 

Reservoir total volume, 10 6 m 3 

12 

12 

12 

n/a 


“Reservoir conditions on testing days are also included. 

I »m 3 denotes I m 3 of gas at 1 bar and 0°C. 

‘Sensitivity analysis range refers to varying gas flow rate while maintaining constant initial plume area. 
‘'Depth at location of lateral CTD transect. 
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a) 2 July 2003 b) 2 July 2003 



-20 -10 0 10 20 -20 -10 0 10 20 
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e ) 23 October 2004 23 October 2004 
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Distance from Diffuser (m) Distance from Diffuser (m) 


Figure 4. Measured (left) plume temperature (°C) and (right) dissolved oxygen (DO) (g/m 3 ) contours 
with linear bubble plume model predictions for diffuser operation with air (2 July 2003) and pure oxygen 
(17 August 2003 and 23 October 2004) in Spring Hollow Reservoir, Virginia. Contours were interpolated 
from CTD profiles collected at locations indicated by small black squares along the bottom of each plot. 


4.3. Sensitivity Analysis 

[ 23 ] A sensitivity analysis was perfonned with the linear 
bubble plume model to determine the effects on plume rise 
height, oxygen transfer efficiency, and induced water flow 
rate, because these parameters are important for design and 
operation of hypolimnetic aeration/oxygenation systems. 
Parameter perturbation was used, in which input variables 
are independently adjusted to determine their individual 
effects on model predictions. The model variables investi¬ 
gated are either difficult to measure or can be controlled 
through system design or operation (Table 3). Currently, a 
and A for the linear plume model are empirical constants 
[Fannel 0 p et al, 1991], and the initial plume water velocity 
is calculated using a densimetric Fr [Wiiest et al, 1992]. 
The initial plume area, gas flow rate, and, to a lesser extent, 
initial bubble size can be controlled through diffuser design 
and operation. The sensitivity of model predictions to 


ambient dissolved nitrogen was also examined for the 
standard case using compressed air (2 July 2003). The 
model assumes that the background dissolved nitrogen 
concentration is equivalent to the saturated value at atmo¬ 
spheric partial pressure and the average hypolimnetic water 
temperature. 

[ 24 ] The DMPR is most influenced by gas flow rate and 
initial bubble radius (Figures 7e and 7f). The gas flow rate is 
directly related to the density of the plume bubble-water 
mixture and subsequently the positive buoyancy and up¬ 
ward momentum of the plume (Table 2). For gas flow rates 
greater than approximately 100 Nm 3 /h, the plume rise 
height is controlled more by the depth of the thermocline 
in SJJR (Figure 5a). As the plume rise height approaches the 
thermocline, further increases in the buoyancy and momen¬ 
tum fluxes can not overcome the strong ambient density 
stratification. The 2003 predictions are more a function of 
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Figure 5. Input boundary conditions for linear bubble plume model validation and sensitivity analysis. 
Data collected from Spring Hollow Reservoir, Virginia, during diffuser operation with compressed air 
(2 July 2003, crosses) and pure oxygen (17 August 2003, solid circles, and 23 October 2004, open circles). 


gas flow rate than those for October 2004 because of 
differing boundary conditions (Figure 5). 

[ 25 ] The plume rise height is moderately sensitive to 
initial bubble radius (Figure 7f). For initial radii less than 
about 1 mm, the DMPRs are nearly independent of bubble 
size. As initial bubble sizes increase, the plumes ascend 
higher and reach a maximum height for a radius of about 
6 mm for each diffuser test. The shapes of the curves can 
be attributed to the dependence of bubble rise velocity 
and gas transfer coefficients on bubble radius [ Wiiest et 
al, 1992]. The DMPR predictions for August are more 
sensitive to initial bubble radius than those for July and 
October (Figure 7f). The plume on 17 August 2003 did 
not have sufficient buoyancy and momentum to rise to the 
thermocline because a relatively low gas flow rate was 
applied (Figure 4c and Table 3). By contrast, the plume 
dynamics for July and October were influenced to a greater 
degree by the thermocline, which in effect damped the 
sensitivity of plume rise height to the initial bubble radius. 

[ 26 ] The predicted DMPR for the linear bubble plume is 
virtually independent of Fr 0 and A over the ranges analyzed 
(Figures 7b and 7c). The plume rise height is moderately 
sensitive to a. Entrainment into the plume is a function of 
plume size, water velocity, and a (Table 1), and entrainment 
of ambient water decelerates the plume. Plume rise is 
influenced to a somewhat greater degree by the initial plume 
area. As the plume area is increased, the buoyancy flux 
decreases because the gas flow rate is constant. Plume rise 
height was unaffected by variations in the ambient dissolved 
nitrogen concentration from 50 to 200% saturation in the 
hypolimnion for the 2 July 2003 diffuser test. Similar to 
the analysis for DMPR, induced water flow rate at the top 
of the plume was found to be insensitive to Fr a and A and 


somewhat more sensitive to a and the initial plume area 
(results not shown). Plume water flow rate is most 
influenced by gas flow rate and initial bubble radius. 
Higher gas flow rates produce greater buoyancy and 
momentum fluxes, which results in greater plume rise 
heights and increased entrainment. 

[ 27 ] Oxygen transfer efficiency (total mass of oxygen 
transferred relative to initial mass of oxygen in bubbles) 
within the hypolimnion was independent of all parameters 
except initial bubble radius, decreasing from nearly 100% to 
around 20% for air and pure oxygen as the initial radii 
increased from approximately 1 mm to 1 cm (Figure 8). 
Oxygen transfer can continue above the DMPR if the 
bubbles are not dissolved. However, secondary plumes are 
not accounted for in the model. Even though undissolved 
bubbles at the top of the plume may continue to transfer 
oxygen during ascent, the oxygen may not be added at the 
desired depth (i.e., below the thermocline). A local mini¬ 
mum with respect to induced water flow rate at the top of 
the plume is predicted for a bubble radius of about 1 mm for 
the conditions in SHR (not shown). This suggests that while 
maximum oxygen transfer efficiency can be achieved with a 
1-mm initial bubble radius, vertical water circulation, and 
hence oxygen distribution in the hypolimnion, will not be 
optimized. As the bubble radius is increased from 1 mm, 
induced water flow rate increases as well for radii up to 
5 mm, but oxygen transfer efficiency decreases rapidly 
within this range (Figure 8). 

5. Comparison of Linear and Circular Bubble 
Plume Models 

[ 28 ] The primary difference between the linear and cir¬ 
cular bubble plume models is the plume geometry. For a 
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Figure 6. In-plume profiles predicted by linear bubble plume model, represented as (a-c) solid lines 
and (d-f) solid lines and symbols. Input data collected from Spring Hollow Reservoir, Virginia, during 
diffuser operation with compressed air (2 July 2003, crosses) and pure oxygen (17 August 2003, solid 
circles, and 23 October 2004, open circles). Measured average in-plume temperature, DO, and plume 
water density represented as symbols (Figures 6a-6c). 


given plume cross-sectional area, the perimeter of a linear 
plume is much greater than for a circular or round plume. 
The initial estimated plume area and perimeter of the linear 
diffuser in SHR for July 2003 was about 50 m 2 and 600 m, 
respectively (Table 3). This corresponds to an equivalent 
radius and perimeter of approximately 4 m and 25 m, 
respectively, for a circular plume. In this case, the perimeter 
of the linear plume is 24 times greater than for the circular 
plume. For both the linear and circular plume models, 
entrainment of ambient water is directly proportional to 
local plume perimeter, local plume water velocity, and a 
(Table 1). The larger perimeter of the linear plume greatly 
increases ambient entrainment, which contributes to the 
negative buoyancy of the plume and causes the plume to 
decelerate more rapidly. This results in a lower plume rise 
height compared with the circular plume. 

[ 29 ] The sensitivity of the linear bubble plume model to 
various parameters was comparable to results for the circu¬ 


lar bubble plume model analysis. Gas flow rate had the 
greatest effect on DMPR predictions by both models, and 
the linear plume model was less sensitive to initial bubble 
radius than the circular plume model was \ Wilest et al., 
1992], The latter may be due to differences between 
temperature boundary conditions, which caused plume 
dynamics in SHR to be more influenced by thermocline 
depth. Initial bubble size greatly affected oxygen transfer 
efficiency for both diffuser geometries, decreasing rapidly 
as the radii increased beyond 1 mm and 3 mm for the linear 
and circular plumes, respectively. Both the linear and 
circular plume model predictions for DMPR were relatively 
insensitive to Fr a and initial plume area (Figure 7 and Wiiest 
et al. [1992]). However, circular model predictions for 
DMPR were more sensitive to a. Overall, the linear model 
was less sensitive to input parameters than the circular 
model, but this insensitivity is probably due to the relatively 
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Figure 7. Effect of linear bubble plume model parameters on depth of maximum plume rise (DMPR). 
Standard values for each parameter are indicated by the vertical dashed lines. Input data collected from 
Spring Hollow Reservoir, Virginia, during diffuser operation with compressed air (2 July 2003, crosses) 
and pure oxygen (17 August 2003, solid circles, and 23 October 2004, open circles). 


homogeneous boundary conditions in the hypolimnion of 
SHR caused by diffuser mixing. 

6. Summary and Conclusions 

[ 30 ] In the current work, the linear bubble plume model 
of McGinnis et at. [2001] was improved, and the updated 
model was evaluated using data collected from a full-scale 
hypolimnetic oxygenation system installed in Spring 
Hollow Reservoir, Virginia. Three diffuser experiments 
were conducted using compressed air and pure oxygen over 
a range of flow rates. Predicted plume rise height, spread¬ 
ing, and constituent profiles were compared with experi¬ 
mental observations. For July 2003 and October 2004, the 
DMPR was simulated well by the model. However, the 
model underestimated the plume rise height for August 

2003, when the gas flow rate was comparatively low. The 
model underpredicted the final plume temperature by ap¬ 
proximately 0.3° and 0.2°C for July 2003 and October 

2004, respectively. The model predicted the plume DO 



Figure 8. Effect of initial bubble radius on oxygen transfer 
efficiency predicted by linear bubble plume model. Input 
data collected from Spring Hollow Reservoir, Virginia, 
during diffuser operation with compressed air (2 July 2003) 
and pure oxygen (17 August 2003 and 23 October 2004). 
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profiles very well for all three diffuser tests, including 
simulating the initial rapid transfer of oxygen immediately 
above the diffuser. 

[ 31 ] A sensitivity analysis was performed to determine 
the effect of various model input parameters. The DMPR 
and induced water flow rate were most influenced by gas 
flow rate and initial bubble radius, moderately sensitive to a 
and the initial plume area, and insensitive to Fr„, A, and 
ambient dissolved nitrogen. Oxygen transfer within the 
hypolimnion was independent of all parameters except 
initial bubble radius, decreasing from nearly 100% to 
around 20% for air and pure oxygen as the initial radii 
increased from approximately 1 mm to 1 cm. 

[ 32 ] The linear bubble plume model for hypolimnetic 
oxygenation has been successfully validated. The model 
can be used to design lake and reservoir oxygenation 
systems and to optimize existing systems to maximize 
oxygen addition. Application of the linear plume model 
requires knowledge of the gas flow rate, initial bubble 
radius, initial plume area, and near-field constituent profiles 
(boundary conditions). Additionally, several empirical 
parameters must be estimated including a and A. Even 
though all of the model inputs may not be known with 
certainty for a given aeration or oxygenation system, the 
model can be used for preliminary design and coarse 
optimization. Also, plume rise height, water flow rate, and 
oxygen transfer efficiency were found to be primarily 
dependent on gas flow rate and initial bubble radius. 

[ 33 ] Operation of bubble plumes for hypolimnetic oxy¬ 
genation usually alters the ambient temperature and DO 
conditions of a water body. Plume dynamics and oxygen 
transfer are strongly related to the near-field water column 
properties, establishing a feedback loop that continually 
changes plume dynamics. This complex plume-lake inter¬ 
action should be accounted for in the design and operation 
of bubble plume diffusers. McGinnis et al. [2001] per¬ 
formed a preliminary coupling of the linear bubble plume 
model with an existing reservoir model, CE-QUAL-W2 
[Cole and Wells, 2003], and obtained encouraging results. 
Efforts are currently under way to further develop the 
coupled model to predict plume performance and the near- 
and far-field reservoir responses. In the absence of near- 
field boundary profiles, far-field or simulated constituent 
profiles can be used to provide a reasonable estimate of 
plume performance. 

Notation 

b plume radius (circular bubble plume), m. 

d bubble diameter, mm. 

C dissolved concentration, mol/m 3 . 

E entrainment factor, m 2 /s. 

F d dissolved species flux, mol/s. 

F c , gaseous species flux, mol/s. 

F s salinity flux, kg/s. 

F t temperature flux, °C m 3 /s. 

Fr Froude number [ ]. 

g gravitational acceleration, m/s 2 . 

H Henry’s constant, mol/m 3 /bar. 

Kl mass transfer coefficient, m/s. 

L plume length, m. 

M water momentum, m 4 /s 2 . 


N number flux of bubbles, 1/s. 

P pressure, bar. 

Q plume flow rate, m 3 /s. 
q actual gas flow rate per unit diffuser length, 
m 2 /h. 

Ri Richardson number [ ]. 

R bubble radius, m. 

5 salinity, g/kg. 

T temperature, °C. 
v velocity, m/s. 

W plume width, m. 
y gaseous concentration, mol/m 3 . 

Z depth, m. 

Greek letters 

a entrainment coefficient [ ]. 

A spreading coefficient [ ]. 
p density, kg/m 3 . 

Subscripts 

3,2 Sauter-mean. 

G Gaussian profile. 

O oxygen. 

N nitrogen. 

T top-hat profile. 
a ambient water. 
b bubble. 

i gas species, oxygen or nitrogen. 
o initial. 

p plume water and gas mixture. 
w plume water. 
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ABSTRACT 

Weiss, R.F., 1974. Carbon dioxide in water and seawater: the solubility of a non-ideal gas. 
Mar. Chem., 2: 203—215. 

New measurements of the solubility of carbon dioxide in water and seawater confirm 
the accuracy of the measurements of Murray and Riley, as opposed to those of Li and Tsui. 
Corrections for non-ideal behavior in the gas phase and for dissociation in distilled water 
are required to calculate solubility coefficients from these sets of data. Equations for the 
solubilities of real gases are presented and discussed. Solubility coefficients for carbon 
dioxide in water and seawater are calculated for the data of Murray and Riley, and are 
fitted to equations in temperature and salinity of the form used previously to fit the solu¬ 
bilities of other gases. 


INTRODUCTION 

Until very recently, direct determinations of the solubility of C0 2 in sea¬ 
water have been limited to a few measurements of Krogh (1904). Most authors 
have preferred to use approximations based on the solubility of C0 2 in aqueous 
sodium chloride solutions. Buch et al. (1932) prepared tables of seawater C0 2 
solubilities from the Bohr (1899) data for sodium chloride solutions, using 
the assumption that the effect of sea salt on the solubility is equal to the ef¬ 
fect of an equal weight of sodium chloride. Lyman (1957) refined this assump¬ 
tion by considering the salting-out of the various constituents of sea salt, sug¬ 
gesting that the Buch et al. values may be about 0.5% too high. 

The lack of direct measurements of C0 2 solubility in seawater has recently 
been alleviated by two independent sets of data: those of Li and Tsui (1971), 
determined by infrared analysis, and those of Murray and Riley (1971), de¬ 
termined gravimetrically. Unfortunately, the agreement between these sets of 
measurements is poor. Whereas the data of Li and Tsui support the values as¬ 
sumed by Buch et al., the solubility data of Murray and Riley are as much as 
3.8% lower than those of Li and Tsui at higher temperatures and salinities. 

The agreement is better at lower temperatures and in distilled water, so that 
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the differences cannot be expressed either as a constant offset or as a con¬ 
stant factor. 

In hopes of resolving these discrepancies, it was decided to make several 
independent measurements of C0 2 solubility. If it could be shown that a few 
accurate measurements covering a wide range of temperature and salinity 
were consistently in close agreement with one of the two sets of data, this 
would provide strong support for the accuracy of that set at all measured 
temperatures and salinities. 

The second objective of this work was the fitting of the best C0 2 solubility 
data to an equation in temperature and salinity of the form used previously 
to fit the solubilities of N 2 , 0 2 , and Ar (Weiss, 1970), He and Ne (Weiss, 
1971a), H 2 (Crozier and Yamamoto, 1974), and Kr (Weiss and Kyser, in prep¬ 
aration). C0 2 is many times more soluble than any of the gases treated previ¬ 
ously, and in the gas phase its departures from the ideal gas approximation 
are large compared to the accuracy with which its solubility can be measured. 
Therefore, before the use of such equations is justified, it is necessary to use 
Henry’s law for real gases and to test its validity over a range of partial pres¬ 
sures of C0 2 . Also, it is necessary to test the validity of the logarithmic Set- 
chbnow salting-out relation (Weiss, 1970,1971b) for C0 2 in seawater. 

Both Li and Tsui (1971), and Murray and Riley (1971), noted that their 
measurements would be affected by the dissociation of dissolved hydrated 
C0 2 to form bicarbonate, and accordingly they acidified their seawater sam¬ 
ples to a sufficiently low pH so that this effect could be neglected. However, 
both studies failed to take this effect into account for their distilled-water 
measurements, which were made without acidification. Although this effect 
cancels in direct comparisons between the two sets of measurements, the 
proper use of Henry’s law and of the Setchbnow relation requires that the 
data be corrected for dissociation. 

All distilled water C0 2 solubility data discussed in the following sections 
have been corrected for dissociation using the following approximation: 

[C0 2 ] = [2C0 2 ] - V K[ [£C0 2 ] (1) 

where [C0 2 ] is the sum of the concentrations of dissolved C0 2 and undis¬ 
sociated hydrated C0 2 , [SC0 2 ] is the sum of the concentrations of all dis¬ 
solved C0 2 species as measured in the solubility experiment, and K{ is the 
first apparent dissociation constant at the temperature of the measurement.* 1 
Solubilities measured in distilled water for p C q — 1 atm are thus reduced 
by ~ 0.18% at 0°C, to ~ 0.46% at 40°C. 

HENRY’S LAW AND REAL GASES 

King (1969, ch.4) presents an excellent and thorough discussion of Henry’s 
law and its application to real gases over a wide range of pressures. By ex- 

* l Values of K\ were calculated according to Hamed and Davis (1943, eq.15): 
log,„k; = —3404.71/T + 14.8435 - 0.032786 T. 
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pressing the activity of the solute in the gas phase by its fugacity, and by 
making the appropriate thermodynamic correction for the expansion of the 
solution by the dissolved gas, the “modified Henry’s law equation” is obtained: 

fi = QiX { exp(PD,/R7) (2) 

where f f is the fugacity of gas i, is the mole fraction of i in solution, v, is 
the partial molal volume of i in solution, R is the gas constant, T is the abso¬ 
lute temperature, P is the total pressure at the liquid-gas interface, and Q { , 
the modified Henry’s law constant for i, is a function only of the temperature 
and the nature of the solvent. 

Conveniently, King uses the CO 2 —water system as an example to illustrate 
the application of eq.2. Using literature data, with special emphasis on the 
work of Wiebe and Gaddy (1939; 1940), King (1969, pp.219—220) demon¬ 
strates that eq.2 holds over the entire C0 2 -pressure range from zero to several 
hundred atmospheres, at a number of different temperatures.* 2 The use of 
the modified form of Henry’s law to calculate C0 2 solubilities over wide 
ranges of total pressures and C0 2 partial pressures, based on accurate solubil¬ 
ity measurements made at pressures near 1 atm, is therefore well justified.* 3 
This is particularly germane to studies dealing with the solubility of natural 
levels of atmospheric C0 2 . 

In treating the solution of atmospheric gases in natural waters, it is con¬ 
venient to express the solubility in terms of the Bunsen coefficient /?, thereby 
avoiding the problem of evaluating the mole fraction of the gas in mixed sol¬ 
vents such as seawater. The Bunsen coefficient is defined here as the volume 
of gas (STP) absorbed per unit volume of the solution, at the temperature of 
the measurement, when the total pressure and the fugacity are both 1 atm.* 4 
This is preferable to the usual definition of P for ideal gases, in which the 
partial pressure is set at 1 atm and the total pressure is not specified, because 
the total pressure is required to define the system (see eq.2). Thus, C,-, the 
volume (STP) of gas i dissolved in a unit volume of solution at the tempera¬ 
ture of the measurement, is given by the relation: 

C,- = Pi fi exp[(l — P)y,/RT] (3) 

where: 

Pi = QJ 1 V} pN exp(—u,/RT) (4) 

and vj is the volume of one mole of the real gas i at STP, p is the density of 
the solution, N is the number of moles in a unit weight of solvent, and P is in 
atmospheres. Since the solubility of most gases (y as 30 cm 3 /mole) is decreased 


* 1 Note that helium solubility measurements by these workers are in excellent agreement 
with recent microgasometric data (Weiss, 1971a). 

* 3 A more detailed discussion of CO, solubility as a function of pressure is given in the 
Appendix. 

* 4 Pressures are given here in atmospheres because of the widespread use of this unit to 
define standard states. Atmospheres may be converted to bars (10 s newtons/m J ), the more 
fundamental unit commonly used in high-pressure work, by multiplying by 0.986923. 



206 


by only ~0.15% per atmosphere increase in total pressure, the small deviations 
from 1 atm total pressure which are encountered in many natural conditions 
may be neglected, thereby reducing the exponential term in eq.3 to unity. 

In the special case of C0 2 , solubilities are generally used in calculations of 
chemical equilibria and are therefore best given in molar (moles/1 of solution) 
or gravimetric (moles/kg of solution) units: 

[CO,] = K 0 f C o 7 exp[(l - P)Uqq JRT] (5) 

J. 4- 

where the constant K 0 equals/3/V 1 in molar units (moles/1 • atm) or (3 /pV in 
gravimetric units (moles/kg • atm). Again, the exponential term may be taken 
as unity at total pressures near 1 atm. 

In the following discussions, gas volumes and fugacity-pressure corrections 
are based on the virial equation of state. The appropriate expansions are carried 
only as far as the second virial coefficient B(T), since the terms containing 
the third virial coefficient were always found to be negligible (~10' s ). Values 
of B(T) in cm 3 /mole, for C0 2 (Sengers et al., 1971) in the range 265—320°K, 
are well represented by a power series: 

B(T) = - 1636.75 + 12.0408 T - 3.27957 • 10' 2 T 2 + 3.16528 • 10~ S T 3 (6) 

Thus, for the pure gas (neglecting the small contribution of water vapor), the 
solubility measurements discussed in the following sections have been corrected 
using the approximations (Guggenheim, 1967, pp.91 and 97): 

V(P, T) = V* (P,T) + B(T) (7) 

and: 

1 n(f/P) = B(T) P/RT (8) 

where V is the volume of one mole of the real gas, and V * is the volume of 
one mole of ideal gas. 

In most natural applications which do not require accuracies greater than 
''■-0.7%, the fugacity in eq.4 may be taken as equal to the partial pressure. 
However, if greater accuracy is desired, the fugacity must be calculated. For 
a nearly pure gas phase, eq.8 will suffice, but for multi-component gas phases, 
such as C0 2 in air, it is necessary to calculate fugacity in the mixture. 

Useful equations for the calculation of fugacities in binary mixtures are 
given by Guggenheim (1967, pp.175—177): 

fl =x 1 Pexp[(fi 11 + 2x\ 5 12 ) P/RT] (9) 

where the subscripts 1 and 2 refer to the two components of the mixture, x is 
the mole fraction, and P n is the second virial coefficient of the pure gas 1. 

The quantity 6 12 is defined by: 

b 12 = 2 + b 22) + 6 12 


( 10 ) 
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where B 12 is the cross virial coefficient for interactions between gas 1 and gas 
2 molecules, and Bn and B 2 2 re ^ er to P ure g 3 * 365 * ^ 2, respectively. In 
addition to the binary mixtures, equations 8 and 9 may also be applied to 
mixtures of C0 2 in air when jcqo 2 « 1* 

Few direct measurements of cross-virial coefficients exist, so that the quan¬ 
tity 6 12 must be evaluated from theoretical considerations. The quantity 
6 co —air has been evaluated from the Lennard-Jones (6-12) potential model 

following the method of Hirschfelder, Curtiss and Bird (1954,^pp.166—170). 
The results in cm 3 /mole for the temperature range 273 to 313°K are well 
represented by the linear relation: 

5 CO,— air = 57.7 — 0.118 T (11) 

The magnitude of this correction, as opposed to assuming 6=0 (Lewis and 
Randall rule), is about 0.2% in the overall solubility calculation, or about one 
third of the total deviation from ideal gas behavior. 

Calculations of C0 2 fugacity using the virial equation of state are valid to 
within 0.1% for total pressures up to ~10 atm. Higher pressures require the 
use of more sophisticated equations of state, such as the Benedict, Webb and 
Rubin (1940; 1942) equation, which may be used to at least 500 atm (see 
Appendix). 

DATA FITTING 

The corrected data of Murray and Riley (1971) were fitted to the same 
equation in temperature and salinity which has been used to fit the solubili¬ 
ties of many other gases (see Introduction). This equation is derived from the 
integrated van’t Hoff equation and the logarithmic Setch^now salinity depen¬ 
dence, and has the form (Weiss, 1970): 

In K 0 = A 1 + A 2 (100/T) + A 3 ln(T/100) + 

^[Bi+BAT/^+BAT/lOO) 2 ] (12) 

where K 0 may be expressed either in moles/1 • atm, referring to a liter of solu¬ 
tion at the temperature of the measurement and an atmosphere fugacity in 
the gas phase, or in moles/kg • atm, referring to a kilogram of solution. The 
A ’s and B's are constants, T is the absolute temperature, and S°ho is the salinity 
in parts per thousand. The method of fitting the data, and the conditions for 
justification of the Setch6now equation and the number of terms used in the 
expansion, are discussed elsewhere (Weiss, 1970). 

The corrected Murray and Riley C0 2 data showed good agreement with 
the equations at all stages of the fitting procedure. The data show a root-mean- 
square deviation from the final fitted equation of 1.4 • 10 " 4 moles/1 • atm in 
K 0 , or about 0.3%. Agreement with the Setch6now relation is shown in Fig.l 
by the random nature of the deviations as a function of salinity for each mea¬ 
sured isotherm. The fitted values of the constants in eq.12 are listed in Table I 
for K 0 in molar and in gravimetric units. Values of K 0 at various temperatures 




Fig.l, Deviations of the Murray and Riley (1971) CO, solubility data (corrected for dis¬ 
sociation and non-ideality) in moles/1 • atm from the fitted curve (eq.12), plotted against 
salinity. Measurements made at constant temperature are connected by lines, and are 
labeled to the nearest °C. The dashed lines show ±1 standard deviation. 


and salinities calculated from these constants are listed in Tables II and III. 

Salting-out constants (the B's in eq.12) were also determined for the solubil¬ 
ity of C0 2 in aqueous NaCl solutions using the data of Bohr (1899), corrected 
for dissociation. With NaCl concentration expressed as weight percent in solu- 

TABLEI 

Constants for the calculation of the solubility of CO, in molar and gravimetric units 
according to eq.12 


Units of K 0 


moles/1 ■ atm moles/kg • atm 


A, -58.0931 

A, 90.5069 

A 3 22.2940 

B, 0.027766 

B, -0.025888 

B 3 0.0050578 


-60.2409 

93.4517 

23.3585 

0.023517 

-0.023656 

0.0047036 
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TABLE II 

The solubility coefficient K 0 (10' 5 moles/1 • atm) 


Salinity ( 0 /«,) 


T(°C) 

0 

10 

20 

30 

34 

35 

36 

38 

40 

-1 

— 

— 

7.273 

6.903 

6.760 

6.724 

6.689 

6.620 

6.551 

0 

7.758 

7.364 

6.990 

6.635 

6.498 

6.465 

6.431 

6.364 

6.298 

1 

7.458 

7.081 

6.723 

6.382 

6.251 

6.219 

6.187 

6.123 

6.060 

2 

7.174 

6.813 

6.469 

6.143 

6.017 

5.986 

5.955 

5.894 

5.833 

3 

6.905 

6.558 

6.229 

5.916 

5.795 

5.766 

5.736 

5.677 

5.619 

4 

6.650 

6.317 

6.001 

5.701 

5.585 

5.557 

5.528 

5.472 

5.416 

5 

6.408 

6.088 

5.785 

5.497 

5.386 

5.358 

5.331 

5.277 

5.223 

6 

Q.178 

5.871 

5.580 

5.303 

5.196 

5.170 

5.144 

5.092 

5.040 

8 

5.751 

5.469 

5.200 

4.945 

4.846 

4.822 

4.797 

4.749 

4.702 

10 

5.366 

5.105 

4.867 

4.621 

4.529 

4.507 

4.485 

4.440 

4.396 

12 

5.017 

4.776 

4.546 

4.327 

4.243 

4.222 

4.201 

4.160 

4.119 

14 

4.700 

4.477 

4.264 

4.062 

3.983 

3.964 

3.945 

3.906 

3.869 

16 

4.412 

4.205 

4.008 

3.820 

3.747 

3.729 

3.712 

3.676 

3.641 

18 

4.149 

3.958 

3.775 

3.600 

3.533 

3.516 

3.499 

3.466 

3.434 

20 

3.910 

3.732 

3.562 

3.400 

3.337 

3.322 

3.306 

3.275 

3.245 

22 

3.691 

3.526 

3.368 

3.217 

3.158 

3.144 

3.130 

3.101 

3.073 

24 

3.491 

3.337 

3.190 

3.050 

2.995 

2.982 

2.968 

2.942 

2.915 

26 

3.307 

3.164 

3.027 

2.897 

2.846 

2.833 

2.821 

2.796 

2.771 

28 

3.138 

3.005 

2.878 

2.756 

2.709 

2.697 

2.685 

2.662 

2.639 

30 

2.983 

2.859 

2.741 

2.627 

2.583 

2.572 

2.561 

2.540 

2.518 

32 

2.840 

2.725 

2.615 

2.509 

2.468 

2.457 

2.447 

2.427 

2.407 

34 

2.708 

2.601 

2.498 

2.400 

2.361 

2.352 

2.342 

2.323 

2.305 

36 

2.587 

2.487 

2.391 

2.299 

2.263 

2.254 

2.246 

2.228 

2.211 

38 

2.474 

2.382 

2.292 

2.207 

2.173 

2.165 

2.157 

2.140 

2.124 

40 

2.370 

2.284 

2.201 

2.121 

2.090 

2.082 

2.074 

2.059 

2.044 


tion, and K 0 in moles/1* atm, the constants are: = —0.68330, B 3 = 0.40911, 
B 3 = —0.064989. Averaged over the temperature range 0—40°C, K 0 is 0.8% 
higher in a 3.6% NaCl solution than K 0 in seawater of 36%, salinity, calculated 
from the salting-out constants in Table I. This difference probably lies within 
the error of the Lyman (1957) prediction that K 0 in the NaCl solution would 
be 0.5% higher. 


EXPERIMENTAL METHOD 

Measurements of C0 2 solubility were carried out by the microgasometric 
technique used previously to measure He and Ne solubilities (Weiss, 1971a). 
High-purity C0 2 (certified > 99.99%) was supplied by Matheson Gas Products 
and gas chromatographic analysis showed < 0.01% air contamination. Because 
of the high solubility of C0 2 compared to He and Ne, the amount of degassed 
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water added to the equilibration chamber was reduced to about one third of 
the total gas volume. Under these conditions, equilibration was extremely 
rapid (first-order rate constant rs 30 sec) so that the 6—10 minutes allowed 
for equilibration were more than adequate. Permeation of C0 2 through the 
indicator drop was rapid, and required that the chamber above the drop be 
filled with pure C0 2 prior to each equilibration. With this procedure, no drift 
in the indicator drop could be detected over a period of 1 h. 

In order to obtain the most direct comparisons with the work of Murray 
and Riley and of Li and Tsui, chemical procedures similar to their’s were fol¬ 
lowed. Distilled water measurements were performed without acidification 
and the corrections for dissociation were made. Surface seawater, collected 
at La Jolla and evaporated to increase its salinity by ~2%o, was passed through 
a 0.45/1 filter, poisoned with 1 mg/liter HgCl 2 , and its salinity determined to 
± 0.004°/oo by an inductive salinometer. Sulfuric acid (~2N) was then added 
to bring the pH to 2.2 and the salinity was adjusted by gravimetric determina¬ 
tion of the amount of H 2 0 added in the acid solution. Following Murray and 

TABLE III 


The solubility coefficient K 0 (10 moles/kg • atm) 


Salinity (7oo) 


T (°C) 

0 

10 

20 

30 

34 

35 

36 

38 

40 

-1 

— 

—- 

7.158 

6.739 

6.579 

6.539 

6.500 

6.422 

6.345 

0 

7.758 

7.305 

6.880 

6.479 

6.325 

6.287 

6.249 

6.175 

6.101 

1 

7.458 

7.024 

6.616 

6.232 

6.085 

6.048 

6.012 

5.941 

5.870 

2 

7.174 

6.758 

6.367 

5.999 

5.857 

5.822 

5.788 

5.719 

5.651 

3 

6.904 

6.506 

6.131 

5.777 

5.642 

5.608 

5.575 

5.509 

5.444 

4 

6.649 

6.267 

5.907 

5.568 

5.438 

5.40.5 

5.374 

5.310 

5.248 

5 

6.407 

6.040 

5.695 

5.369 

5.244 

5.213 

5.182 

5.122 

5.062 

6 

6.177 

5.825 

5.493 

5.180 

5.060 

5.031 

5.001 

4.943 

4.885 

8 

5.752 

5.427 

5.120 

4.831 

4.720 

4.693 

4.666 

4.612 

4.558 

10 

5.367 

5.067 

4.784 

4.516 

4.413 

4.388 

4.363 

4.313 

4.263 

12 

5.019 

4.741 

4.479 

4.231 

4.136 

4.112 

4.089 

4.042 

3.997 

14 

4.703 

4.446 

4.202 

3.972 

3.884 

3.862 

3.840 

3.797 

3.755 

16 

4.416 

4.177 

3.951 

3.738 

3.655 

3.635 

3.615 

3.575 

3.536 

18 

4.155 

3.933 

3.723 

3.524 

3.448 

3.429 

3.410 

3.373 

3.336 

20 

3.916 

3.710 

3.515 

3.330 

3.258 

3.241 

3.223 

3.189 

3.154 

22 

3.699 

3.507 

3.325 

3.152 

3.086 

3.069 

3.053 

3.021 

2.989 

24 

3.499 

3.321 

3.151 

2.990 

2.928 

2.912 

2.897 

2.867 

2.837 

26 

3.317 

3.150 

2.992 

2.841 

2.783 

2.769 

2.755 

2.727 

2.699 

28 

3.149 

2.994 

2.846 

2.705 

2.651 

2,638 

2.624 

2.598 

2.572 

30 

2.995 

2.850 

2.712 

2.580 

2.530 

2.517 

2.505 

2.480 

2.455 

32 

2.854 

2.718 

2.589 

2.466 

2.418 

2.406 

2.395 

2.372 

2.349 

34 

2.723 

2.596 

2.476 

2.360 

2.316 

2.305 

2.294 

2.272 

2.250 

36 

2.603 

2.484 

2.371 

2.263 

2.221 

2.211 

2.201 

2.180 

2.160 

38 

2.492 

2.381 

2.275 

2.174 

2.134 

2.125 

2.115 

2.096 

2.077 

40 

2.389 

2.285 

2.186 

2.091 

2.054 

2.045 

2.036 

2.018 

2.000 
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Riley (1971), the effect of the added H 2 S0 4 molecules on the salinity was 
neglected as being insignificant with respect to altering the solubility of C0 2 . 
The water was degassed using the same vacuum extraction method as in the 
previous work (Weiss, 1971a). 

The results were calculated in the same manner as previous microgasometric 
measurements, except that the corrections for non-ideal behavior discussed in 
the foregoing section were applied. The correction for the expansion of the 
aqueous phase during equilibration was made using a partial molal volume of 
32.3 cm 3 per mole of C0 2 (see Appendix). The overall accuracy of these C0 2 
solubility measurements is estimated as ±0.2%. 

RESULTS AND CONCLUSIONS 

The microgasometric experimental values for K 0 in moles/1 • atm at several 
temperatures and salinities, corrected for dissociation and non-ideality, are 
given in Table IV. Deviations of the microgasometric results from the fitted 
Murray and Riley data, plotted in Fig.2, show excellent agreement without 
systematic deviations. At each of the three different temperature and salinity 
conditions, the range of microgasometric points brackets the fitted curve. The 
average deviation of all the measured points is 1,5 • 10 _s moles/1 • atm, and the 
root-mean-square deviation is 1.2 • 10 -4 moles/1 • atm. Curves for the deviations 
of the data of Li and Tsui, taken from their equations (I) and (2a) after cor¬ 
rection for dissociation and non-ideality, are also plotted in Fig.2. Although 
their fresh-water measurements are in reasonably good agreement, the seawater 
data of Li and Tsui show large deviations from the fitted Murray and Riley 


TABLE IV 


Experimental C0 2 solubilities: microgasometric determinations of K 0 in moles/1-atm 


Salinity (%o) 

f(°C) 

K 0 -10 s 

0.0 

20.61 

3.848 

0.0 

20.63 

3.843 

0.0 

20.60 

3.867 

0.0 

20.59 

3.834 

0.0 

20.61 

3.840 

35.330 

6.59 

5.071 

36.330 

6.60 

5.063 

35.330 

6.60 

5.068 

35.330 

6.59 

5.041 

35.330 

20.62 

3.245 

35.330 

20.63 

3.272 

35.330 

20.63 

3.245 

35.330 

20.63 

3.255 

35.330 

20.63 

3.254 

35.330 

20.63 

3.252 
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Fig. 2. Deviations of various CO, distilled water and seawater solubility data in moles/l .* atm 
from the fit to the corrected data (see text) of Murray and Riley (1971), plotted against 
temperature. The curves for Li and Tsui (1971) are taken from their equations (1) and 
(2a), and, like the data of Krogh (1904), have been corrected for non-ideality and dissocia¬ 
tion: the dashed portions of these curves represent the extrapolation of their equations 
beyond the range of their measurements. The magnitude of the error caused by dissociation 
for measurements made in distilled water at pCO, = 1 atm is shown by the line labeled HCO, . 


data. Fig.2 also shows the corrected data of Krogh (1904), which agree well 
with Murray and Riley at the higher temperatures. 

The microgasometric results provide independent confirmation of the accu¬ 
racy of the CO 2 solubility measurements of Murray and Riley (1971) and 
show the measurements of Li and Tsui (1971) to be in error by as much as 
~4% at the higher temperatures and salinities. The data of Murray and Riley, 
corrected for the effects of non-ideal gas behavior and for dissociation in the 
distilled water measurements, are well represented by the Setch6now and 
integrated van’t Hoff equations used previously to fit the solubilities of several 
other comparatively ideal gases. Thus, C0 2 solubilities calculated from eq.12, 
using the constants in Table I, are believed to provide the most accurate values 
in the literature, with an overall accuracy estimated at ±1 • 10" 4 moles/1* atm 
or about 0.2%. The solubility of C0 2 obeys the modified form of Henry’s law 
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(eq.2) and it is therefore necessary to take account of the total pressure, as 
well as to calculate the fugacity of C0 2 in the gas phase, if full use is to be 
made of the accuracy of these solubility data. 
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APPENDIX 

Carbon dioxide solubility at high pressures 

At pressures greater than ~10 atm, CO, fugacities calculated from the virial equation 
of state, even when carried to the third virial coefficient, show large deviations from the 
literature values (Deming and Deming, 1939). Far better results are obtained with the 
Benedict, Webb and Rubin equation of state for pure substances (Benedict et al., 1940) 
and for mixtures (Benedict et al., 1942), which fit the observed p—V—T data up to densi¬ 
ties of twice the critical density. In the following discussion, Benedict, Webb and Rubin 
equation constants for pure CO, are taken from the Bishnoi and Robinson (1971 )* 5 
volume-dependent fit to the compressibility data of Reamer et al. (1944). Because the 
Benedict, Webb and Rubin equation is an explicit function of V and T, the solution for 
specific values of P is obtained by the method of successive approximations. 

Fig.3 shows the distilled water COj solubility data of Wiebe and Gaddy (1940) for the 
temperature range 12—40°C and the pressure range 25—500 atm, plotted as In (f/x) against 
the total pressure P. Fugacities were calculated assuming a pure CO, phase at a pressure of 
P minus the vapor pressure of water. Negligible error was introduced by this approximation 
because of the small fraction of water vapor in the CO,-rich phase. The straight lines plotted 
at each of the six experimental temperatures were fitted to the Wiebe and Gaddy data by 
the method of least-squares. 

According to eq.2, if v is independent of P, then In (f/x) plotted as a function of P at 
constant T should give a straight line with intercept In Q and slope (C/RT). This relationship 
is strongly supported by the experimental data as shown in Fig.3. The root-mean-square 
deviation of the Wiebe and Gaddy measurements about the six straight lines is ~0.8%, and 
there is no systematic indication of departure from linearity. The slopes of these lines cor¬ 
respond to remarkably constant values of v, with a mean of 32.3 cm 3 /mole and a standard 
deviation of 0.5 cm 3 /mole. Thus, the data show no significant variation in v, either as a func¬ 
tion of T over the range 12—40°C, or as a function of P up to 500 atm. 

Results of the fit to the corrected Murray and Riley distilled water data (eq.12), mea¬ 
sured at 1 atm, are also plotted in Fig.3. The overall agreement with the Wiebe and Gaddy 
data is good, although there may be a systematic difference at 12°C (the Wiebe and Gaddy 
data also show the greatest scatter at this temperature). On the average, the solubilities at 
1 atm obtained from the linear fits to the Wiebe and Gaddy data are 0.4% lower in 0 (or 
higher in Q), and show a 2% root-mean-square difference. 


* s The constants, converted from English to metric units (atm, 1, mole, °K): A 0 = 1.9521940, 
B„ = 0.033065171, C 9 = 170449.55, a ■= 0.25264437, b = 6.6041298 • 10' 3 , c = 19592.012, 
a = 4.7117009-10' 5 , y = 4.3414415• 10‘ 3 , R = 0.08205601. 
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l-—-1-1-1-!----L.j 

o 100 200 300 400 500 


P(atm) 

Fig.3. High-pressure distilled water C0 2 solubility data of Wiebe and Gaddy (1940) plotted 
as the logarithm of the fugacity to mole fraction ratio, against the total pressure. Solubili¬ 
ties at 1 atm calculated from eq.12 are plotted for the temperatures of the Wiebe and 
Gaddy measurements. 


In summary, the high-pressure solubility data show that the modified form of Henry’s law 
(eqs.2, 3, or 5) is valid for GO, in water over the entire pressure range of 0—500 atm, Which 
corresponds to a dissolved CO, concentration range of 0—1.7 molar. Assuming that v for 
CO, is the same in seawater as in distilled water, as was shown by Enns et al. (1965) for O,, 
the solubility of CO, in seawater at high pressures may also be calculated from these equations. 
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Ellipsoidal Fluid Particles I IX. Fluid Dynamics 
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systems from the resulting correlation. Cases where wall effects were too sig¬ 
nificant were also eliminated from the data treated. The criteria which the data 
had to meet were then Eqs, (9-33) or (9-34) and 




The indices in the original Johnson and Braida correlation and the point of 
intersection between the two linear regions were adjusted to improve the agree¬ 
ment with all the data meeting the above criteria. The resulting correlation 
(G12) is: 

./ = 0.94W 0 ' 5? (2 < H < 59.31 (7-51 
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Ellipsoidal Fluid Particles II. Fluid Dynamics 
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